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ABSTRACT 


The development of appropriate modeling and adjustment 
procedures for the estimation of harmonic coefficients of the 
geopotential, from surface gravity data was studied, in order to 
provide an optimum way of utilizing the terrestrial gravity information 
in combination solutions currently developed at NASA/Goddard Space 
Flight Center, for use in the TOPEX/POSEIDON mission. 

The mathematical modeling was based on the fundamental boundary 
condition of the linearized Molodensky boundary value problem. 
Atmospheric and ellipsoidal corrections were applied to the surface 
anomalies. The l’xl* mean free-air anomalies of the Ohio State 
University’s June 1986 global anomaly field were used. Low degree 
(<50) potential coefficient sets were estimated through a rigorous least 
squares adjustment. The high frequency content of the anomalies was 
removed prior to the adjustment using the OSU86F high degree 
expansion. 

Terrestrial gravity solutions were found to be in good agreement 
with the satellite ones over areas which are well surveyed 
(gravimetrically), such as North America or Australia. However, 
systematic differences between the terrestrial only models and GEMT1, 
over extended regions in Africa, the Soviet Union and China, were 
found. In Africa, gravity anomaly differences on the order of 20 mgals 
and undulation differences on the order of 15 meters, over regions 
extending 2000 km in diameter, occur. Comparison of the GEMT1 
implied undulations with 32 well distributed Doppler derived 
undulations gave an RMS difference of 2.6 m, while corresponding 
comparison with undulations implied by the terrestrial solution gave 
RMS difference on the order of 15 m, which implies that the terrestrial 
data in that region are substantially in error. 
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CHAPTER I 


INTRODUCTION 

Accurate knowledge of the external gravity field of the earth is a 
prerequisite for various geodetic and geophysical investigations and 
applications. It is essential for accurate determination of the orbits of 
artificial satellites) since the earth’s gravitational perturbations are 
(for common satellite altitudes) the most significant ones experienced 
by the satellites. Accurate gravity field models are also important for 
geophysical investigations of the earth’s interior, and for 
oceanographic studies related to ocean dynamics. 

The estimation of the external gravity field of the earth, on a 
global basis, has been investigated in the past from the theoretical as 
well as the computational point of view. Improvement in our knowledge 
of the gravity field has become a continuous effort, primarily due to 
the following reasons: 

- the quality, quantity and distribution of the necessary data for 
accurate global gravity modeling, is improving with time as well 
as with the development and use of more accurate observational 
systems as the modern, highly accurate and stable laser tracking 
systems and the altimeter satellites. 

- the availability of more accurate data requires more rigorous and 
complete mathematical modeling, since the effects of certain 
approximations made in the past become non- negligible in the 
presence of more accurate and complete data. 
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- some of the approximations which were applied in order to make 
the estimation of global gravity models a computationally 
manageable problem, become unnecessary in view of the 
computational resources available nowadays, especially after the 
development of vector processors. Estimation techniques which 
have found very limited or no application in the paBt due to 
their computational load, can now be implemented and their 
possible merits can be investigated. 

Currently, the estimation of global gravity field models is based on 
the combination of information coming from the following main sources: 

(a) Artificial satellite observations, such as laser ranging, radar 
and optical observations, 

(b) Terrestrial gravity measurements, primarily gravity anomalies, 

(c) Satellite altimeter observations. 

Satellite observations provide a global sampling of the gravity 
field, but lack detailed information due to the attenuation of the field 
at satellite altitudes (LAGEOS, for example, is practically insensitive to 
harmonics of degree larger than eight). For this reason, satellite 
observations can support accurate determination of the low frequency 
part of the spectrum of the gravity field, but are inadequate for 
estimation of the medium and high frequency part of it. 

The opposite is true for terrestrial measurements, which provide 
detailed sampling of the gravity field, but lack global coverage due to 
the existing data gaps. 

The altimeter observations constitute a special case of satellite 
observations which provide a quite detailed and accurate sampling of 
the field over the oceanic areas but suffer from the large gaps 
generated by the land discontinuities. Therefore, their character and 
information content correspond, as far as the estimation of the gravity 
field is concerned, more to the terrestrial measurements than to the 



3 


satellite ones. 

Consequently, global gravity models developed for "general 
purpose", i.e. both for terrestrial use and for satellite orbit 
computations, need to utilize both satellite and terrestrial data sources 
in a complementary way, so that an optimum determination of the 
gravity field over a wider band of its spectrum could be achieved. At 
this point "special purpose" gravity models such as PGS-S4 [Lerch et 
al, 1982b] "tailored" to fit the orbits of particular satellites (in the 
above case SEASAT) or the Rapp 1977 model [Rapp, 1977], computed 
from terrestrial data alone will not be considered. 

1.1 Previous Investigations 

The necessity to determine optimum ways of combining the satellite 
information with the terrestrial data, for the estimation of global 
gravity models, has been recognized early in the past. Investigations 
by Kaula (1966b) and Rapp (1967) resulted in the formulation of two 
main estimation techniques. The basic difference between the two 
methods is the way that the terrestrial gravity data are being treated. 
In Kaula’B approach the orthogonality relations are used in order to 
determine potential coefficients from gravity anomalies. In Rapp’s 
approach the problem is essentially inverted and a least squares 
adjustment is used to fit a truncated set of potential coefficients to 
the given gravity anomalies. 

Rapp (1969) haB compared both analytically and numerically the two 
estimation techniques and has investigated the conditions under which 
they yield similar results. 

The major disadvantage of Kaula’s approach as opposed to Rapp’s 
is that the gravity anomalies need to be reduced to a spherical surface 
so that the orthogonality relations can be implemented. On the other 
hand, in Rapp’s method a truncated set of coefficients is used to fit 
the data; consequently, aliasing effects will occur (in general) if the 
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data contain frequencies beyond the maximum frequency implied by the 
truncation [Desrochers, 1971]. Attempts to overcome this problem by 
extending the expansion to higher degrees result in large and 
computationally unmanageable normal matrices. 

A detailed comparison of these two techniques will be presented in 
Chapter IV. 

An alternative approach to the problem was proposed and tested 
by Colombo (1981). His method uses again the orthogonalities but now 
a linear estimator, optimal in the sense that it minimizes the sum of 
squares of the coefficient errors, is used in order to provide estimates 
of the potential coefficients. 

These main estimation methods (with some modifications) have been 
implemented in the past by several investigators to provide combination 
solutions for potential coefficient expansions to various degrees. In 
Table 1 some of these solutions which are currently used in several 
studies are presented. 


Table 1. Some Currently Used Global Geopotential Models 


Name 

Author 

Date 

N 

l 'ma x 

Rapp78 

Rapp 


180 

GEM10B 

Lerch et al 


36 

GEM10C 

Lerch et al 


180 

Rapp81 

Rapp 

1981 


Hajela84 

Hajela 

m 

250 

GRIM3-L1 

Reigber 


36 

GPM2 

Wenzel 


200 

0SU86C/D 

Rapp and Cruz 


250 

0SU86E/F 

Rapp and Cruz 

1986 

360 
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In view of the wide applications that such gravity models have 
found and the increasing demand for accurate high degree potential 
coefficient expansions, it has become quite critical to investigate and 
reconsider both the modeling and estimation part of the combination 
solutions. 

The implementation of refined models to account for the effects of 
the earth's ellipticity, the atmosphere and the terrain has been 
investigated by Rapp (1977, 1984, 1986), and special emphasis was 
placed on the problem of the continuation of surface gravity data to a 
spherical bounding surface. Cruz (1985, 1986) has also studied the 

effects of the earth's ellipticity and of the downward continuation and 
his proposed approach has been implemented in recent high degree 
expansions at The Ohio State University [Rapp and Cruz, 1986a, b]. 

Although (as it is evident from the above) the various aspects of 
the determination of gravity models have been extensively studied in 
the past, it is not reasonable to believe that all the problems have 
found satisfactory answers and that there is no possibility of further 
improvement. Rather, the previous investigations and experiences have 
raised new questions that warrant careful study. 


1.2 Objective and Organization of the Present Study 

Currently, an effort to provide a highly accurate global gravity 
model, to be used for the TOPEX/POSEIDON mission [JPL, 1985], is 
undertaken by three main research centers: 

(1) NASA/Goddard Space Flight Center (GSFC) 

(2) University of Texas at Austin 

(3) The Ohio State University. 
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The task here is to investigate optimum modeling and estimation 
procedures for the processing and the analysis of the terrestrial 
gravity information. This information, in the form of a system of 
normal equations referring to a low degree (<50) harmonic expansion of 
the geopotential, will then be combined at NASA/GSPC with the 
corresponding normal equations obtained from the analysis of satellite 
observations, in order to provide an optimum combination solution for 
the potential coefficients. 

The potential coefficient estimates obtained from the terrestrial 
gravity data alone, although they do not constitute the final product of 
this work, are obviously of great value for comparison purposes, in 
order to investigate possible problems related to the surface gravity 
data and to indicate necessary improvements in the modeling or 
estimation techniques used. 

The investigation which was conducted in order to accomplish the 
purpose of this study, can be divided into four main parts as follows: 

(a) Theoretical investigation of the appropriate mathematical models 
that should be used in order to relate the surface mean free-air 
anomalies to potential coefficients. The intention here is to establish 
rigorous mathematical models and avoid any kind of unnecessary 
approximations. These aspects are presented in Chapter II in an 
explicit and detailed form, following a step-by-step discussion which, 
starting from the formulation of Molodensky’s boundary value problem 
(which constitutes the underlying theoretical background of the 
problem at hand), leads to the final analytical formulation used in the 
estimation procedure. 

(b) Description of the available data. The plan in this study is to use 
the most up to date gravity and elevation information that is available. 
In Chapter III the description of the main data sets that were examined 
and/or used in this study is given, focusing on some of their 
characteristics that should be taken into account for the appropriate 
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use of these data and the accurate interpretation of the results. 

(c) The estimation procedure. In this study the estimation of potential 
coefficients from surface gravity is to be made through a least squares 
adjustment procedure as it was proposed by Rapp (1967). In Chapter 
IV the detailed formulation of this procedure is given. In addition its 
properties are examined and a comparison of this technique to the 
quadrature formulas is presented. 

(d) The numerical analysis. In Chapter V the results of a number of 
experiments performed, are presented and discussed. These include 
both the results of experiments performed with simulated data, in order 
to examine several aspects of the procedure, and the results obtained 
from actual data. The latter are compared with existing geopotential 
models in order to validate the solutions and investigate the agreement 
of the gravity models implied by terrestrial data alone, with the 
satellite only or the combined geopotential models. 

Finally, in the last chapter a summary of this work, a presentation 
of the conclusions drawn from it, and recommendations for further 
investigation are given. 


CHAPTER II 


MODELING OF THE TERRESTRIAL GRAVITY DATA 


2.1 Molodensky’s Boundary Value Problem 

The determination of the geometric shape and the external gravity 
potential of the earth from observations carried out on its unknown 
surface, is a non-linear free boundary value problem (bvp), in the 
sense that the geometric shape of the boundary surface is also 
unknown. This geodetic bvp or Molodensky’s problem can be 
formulated as follows: 

Given at all points of the physical surface of the earth (S) the 
gravity potential W and the gravity vector find the shape of the 
surface (S), as well as, the potential W outside (S), so that W consists 
of a regular at infinity harmonic part V, satisfying Laplace’s equation 
* 2 V = 0 and a known centrifugal part, and, satisfies known boundary 
conditions on (S). 

This type of problem can be considerably simplified after 
application of a linearization procedure. Linearization, in the context 
of the modern approach of Molodensky, consists of the introduction of 
a known surface (E), the telluroid, which is "sufficiently close" to the 
physical surface of the earth (S), and the introduction of a normal 
gravity potential U as an analytical approximation of the actual gravity 
potential W. A point P on (S) (Figure 1) will be mapped to a point Q 
on (E) through a one-to-one correspondence. The particular rules that 
can be used to define this mapping will be discussed in Section 2.1.2. 
The important issue here is that (E) is defined in such a way that the 
value of a function on (S) can be obtained from the corresponding 
value of this function on (E) by applying a linear correction . If 
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= Q?, then a Taylor series expansion in terms of ^ around (E) can be 
performed, where all second and higher order terms in ^ can be 
neglected without any significant loss of accuracy (see also the 
discussion at the end of Section 2.1.1). 

After linearization the geodetic bvp becomes a linear non-free bvp 
for the disturbing potential T=W-U, which satisfies Laplace's equation 
outside (E) with known boundary conditions on (E). Once this problem 
is solved, the deviation ^ of (E) from (S) can be determined and 
consequently both the external gravity potential W = U + T and the 
shape of the earth (S) can be defined. 

The approach outlined above is in contrast to the traditional 
approach of Stokes where gravity observations on the earth’s surface 
are first reduced to the geoid, which in turn is approximated by a 
rotational ellipsoid (spheroid) or a sphere. 


2.1.1 The Fundamental Boundary Condition 

Krarup (1973) has rigorously linearized the free non-linear geodetic 
bvp. We are presenting in the following an outline of his approach 
following Moritz’s (1980) notations. 

The basic geometry associated with our problem is described in 
Figure 1. 

We proceed with some basic definitions: 

W = gravity potential of the earth 

U = normal gravity potential defined through an ellipsoid of 

revolution whose surface is an equipotential surface of its 
gravity field (Somigliana-Pizzetti normal field) [Heiskanen 
and Moritz, 1967, Sec. 2-7]. 




11 


?p = + M a? 

where 


(2.9) 


Mn — 


3 2 U 

3 2 U 

a 2 u 

3x a 

3x3y 

3x3z 

3 2 U 

a 2 u 

a 2 u 

3y3x 

ay 2 

3y3z 

a 2 u 

3 2 U 

a 2 u 

dzdX 

3z3y 

3z 2 


( 2 . 10 ) 


is the (second-order) normal gravity gradient tensor (Marussi tensor). 
From (2.2) and (2.8) we have 

Up = U 0 + •? (2.11) 


Similarly, 

T P = T q + (gradT)„-? (2.12) 

However, since gradT*^ is 0(f 2 ) it can be neglected in the context of 
linearization, so that 

Tp = T q (2.13) 

Omissions as the above will be uniformly applied in the following 
derivations. In other words, no distinction will be made between (S) 
and (2) as far as quantities of the anomalous field such as T, Ag, 
gradT, ^ etc. are concerned [Moritz, 1980, pp. 339-341]. Hence, 


AW = Wp - U Q 

= Tp + Up - U q 
- Tq + Uq + ^q'{ ~ Uq Or 

AW = Tq + (2.14) 


and 


Ag = tp _ ?Q 

- tp ~ fp + Mq? 

= [grad(W-U)] P + M q ? 


or 
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Ag = (gradT)p + Mg? 

and since (gradT)p * (gradT)g one finally has 

Ag = (gradT)q + Mq? (2.15) 

Assuming that det(M 0 ) * 0, so that M Q _l exists one has from (2.15) 

? = M 0 _l [Ag - (gradT)g] (2.16) 

so that (2.14) can be written as 

T q + Sq- ( gradT)q = AW + ifcq-Ag (2.17) 

where 

ift Q = -Mq- l f Q (2.18) 

defines a vector tangent to the isozenithal lines [Moritz, 1980, p. 345] 
of the normal gravity field. Equation (2.17) is the fundamental 
boundary condition of the linearized Molodensky bvp. 

Introducing the orthogonal curvilinear coordinates />, ♦, X where 4» 
is the normal latitude, X the normal longitude (see also [Moritz, 1980, p. 
338]) and p is measured along the isozenithal, we can define a set of 
generalized coordinates as 


7i 


- — T cos+cosX 



p 

72 

— 

- "“*r cos+sinX 



P 

73 


- sin* 



P J 


(2.19) 


With respect to these coordinates, equation (2.17) takes the form 
[Moritz, 1980, pp. 342-348] 


AW 


if = -[ill + — [i^[ T a - — M 

1*tJq 7q 1 «tJ q Q 7q l«rJ q 


( 2 . 20 ) 
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where 

y is the magnitude of the normal gravity vector y 
t is the arc length along the isozenithal line 

Ag' is the magnitude of the projection of Ag along the isozenithal 
line passing through Q. 

According to the above, the procedure for the determination of (S) 
and W can be idealistically outlined as follows: 

(a) At every point P on (S) gravimetric and astronomical observations 
provide the magnitude and direction of the gravity vector ^p. On 
the other hand, leveling combined with the gravimetric 
observations provide the gravity potential W up to an unknown 
constant that can be determined indirectly by distance 
measurement (s) [Heiskanen and Moritz, 1967, Sec. 2-19, 2-20]. 

(b) At every corresponding point Q on the telluroid (E) the magnitude 
and direction of the normal gravity vector can be computed 
since both the surface (E) and the normal gravity potential are 
known. 

(c) The direction of the isozenithal line passing through Q can be 
determined, once the normal potential is defined. Hence the 
magnitude of the isozenithal projection Ag' , of Ag, can be 
determined. 

(d) Once the mapping P -* Q, that defines the telluroid (E) is 
established, the quantity AW in (2.20) will be defined. Hence, (2.20) 
takes the form 

c i(|^) Q + C 2 T q = C 3 (2.21) 

with known constant coefficients 
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C, = -l 



C 3 = Ag' + C a AW 


( 2 . 22 ) 


(e) Disregarding the fact that the telluroid can be below the surface 
of the earth [Moritz, 1980, p. 355], and assuming that the effect of 
the mass of the atmosphere has been taken into account [ibid, p. 
442], we can consider the disturbing potential T to be harmonic 
outside (2), i.e. 

v 2 T = 0, outside (2). (2.23) 

Consequently, the problem here is to solve (2.23) for T, subject to 
the boundary condition (2.21). This corresponds to a third bvp of 
potential theory (Robin’s problem - see [Lebedev et al., 1979, p. 27]) 
although it is considerably more complicated because it involves oblique 
derivatives (the isozenithal lines are not normal to the telluroid (2)). 
In any case once this bvp is solved for T, then ^ (and thus the 
geometric shape of S) can be obtained from equation (2.16). 


The steps outlined above are meant to describe the principles 
behind the formulation of Molodensky’s bvp. From the practical point 
of view, the prerequisites of Molodensky’s problem, especially the 
continuous coverage of the whole earth’s surface with gravity 
measurements raise serious questions about the relevancy of the 
problem itself from the geodetic point of view [Moritz, 1980]. 


Although equation (2.20) is theoretically rigorous, it is not 
convenient for practical purposes, due to the presence of the 
isozenithal in both the definition of the directional derivatives 3/3r, 
and of Ag'. However, since the curvature of the normal plumb line is 
very small and regular [Heiskanen and Moritz, 1967, p. 196; Moritz, 
1983, p. 7], the normal plumb line can be considered to coincide with 
the straight ellipsoidal normal, in which case the isozenithal will also 
coincide with the straight ellipsoidal normal [Moritz, 1980, pp. 345-346], 



Accordingly, if 
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isozenithal * normal plumb line = straight ellipsoidal normal (2.24) 


equation (2.20) becomes 


Ag' = 





(2.25) 


where Ag' is hereon to be considered the magnitude of the ellipsoidal 
normal projection of Ag and h denotes the arc length along the 
ellipsoidal normal. 


Consider now, <1$, <t H } P the unit vectors on the three axes of 

the local astronomic frame (east, north, zenith) at P, and nt\, <t$, <l h } p 
the corresponding unit vectors for the local normal frame (Figure 2); 
then 

(£„) P = (f £* + t| £x + cos0 t h ) P (2.26) 

where 9 is the total deflection of the vertical and £, t] are its 
latitudinal and longitudinal components respectively [Heiskanen and 
Moritz, 1967, fig. 2-13]. 



Figure 2. The Local Astronomic and Geodetic Frames. 


The gradient of the gravity potential with respect to the local 


normal frame is 


(gradW)p = (j^ £* + 


*W fW £ 1 

Ncos^ax &X #h h J | 


(2.27) 


where M and N are the meridional and prime vertical radii of curvature 
respectively. Accordingly, 


Ag' = -<f P - ?„)•£„„ 

Ag' = -(gradW) P -4 Q ~ Ifql 


(2.28) 


If we assume that point Q lies on the straight ellipsoidal normal 
through P then 


and (2.28) yields 


-Irql = Ag' + (^) 


(2.29) 


(2.30) 


l^pl = -(gradW)p-tu 


= (m«* + Ncos<MX + «h + 71 + cos0 *h)l 


,4 i f - aw aw . awi 

1 ^ M«+ 71 Ncos<MA cos ahi 


(2.31) 


Adding (2.30) and (2.31) by parts one finally gets 


•tpl - I^qI = Ag' + [(1-cosO) ^ " * m " V 


aW 1 

Ncos+axJ p 


(2.32) 


Ag = Ag' + c| 


with 


aw aw 

= l (1 -° os0) ih " ( am ~ 71 


Ma$ 1 Ncos^axJ p 


(2.33) 


(2.34) 


and Ag = Ijfpl - I^qI being the usual gravity anomaly as obtained from 
surface observations of the magnitude of the gravity acceleration. 



According to the above, equation (2.25) becomes 
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Ag 





AW + cp 
Q 


(2.35) 


which is the form of the boundary condition of the linearized version 
of the geodetic bvp, that will be used in the sequel. Jekeli (1981), 
using the same assumptions as above, arrived at the form (2.35) of the 
fundamental boundary condition using a Taylor series expansion in f 
and neglecting second and higher order terms. 


It should be noted here that e p depends on the mapping used to 
define the telluroid, since to derive (2.34) one had to assume (2.29) 
which restricted Q to lie on the straight ellipsoidal normal through P 
and thus has removed two of the three degrees of freedom for the 
determination of the location of Q. 


Concluding this section some remarks should be made about the 
validity of the linearized version of the geodetic bvp. As it is shown 
by Jekeli [1981, eq. 4.6, 4.7] linearization causes an error on the order 
of 


_1 

70 



(2.36) 


for the height anomaly, and on the order of 


< 2 - 37) 

for the gravity anomaly. Using nominal values one has an error on 
the order of 3 mm in the height anomaly and 1.5 /jgals in the gravity 
anomaly. Hence, in view of the present accuracies, linearization is 
completely justifiable. 
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2.1.2 Definition of the Telluroid 

In the previous discussion the mapping that defines the telluroid 
was left unspecified. The general requirement was that this surface 
should be sufficiently close to the earth’s surface, so that linearization 
could be performed. However, as was mentioned, the form (2.35) of the 
boundary condition that will be used, implicitly restricts part of the 
definition of the telluroid point Q. 

The introduction of different mappings to define the telluroid stems 
from the fact that for certain mappings the form of the boundary 
condition becomes simpler. For the Marussi mapping [Moritz, 1980, p. 
338] defined by 

U Q = W P and (2.38) 

and the gravimetric mapping [ibid, p. 339] 

ip =?Q 

the potential anomaly AW and the gravity 
respectively zero, making the right side 
(2.17) simpler. The mapping that will be 
follows: 

If 

C P = W 0 - W P (2.40) 

is the geopotential number at P [Heiskanen and Moritz, 1967, p. 56], 
then the telluroid point Q that corresponds to P lies on the ellipsoidal 
normal through P and is such that 

U Q = U 0 - Cp (2.41) 

where U 0 is the potential on the surface of the equipotential reference 
ellipsoid. 


defined by 

(2.39) 

anomaly vector Ag become 
of the boundary condition 
adopted here is defined as 


According to this definition the potential anomaly becomes 
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AW = W P - U Q 

= W 0 - C P - U 0 + C p or 

AW = W 0 - U 0 = const. (2.42) 

i.e. represents the constant difference between the potential on the 
geoid minus the potential on the surface of the reference ellipsoid. 
Even if one assumes that W 0 = U 0 , this mapping is different from the 
Marussi mapping because here 



Rather, one can say (neglecting the curvature of the normal plumb 
lines) that in this mapping 



The form (2.35) of the boundary condition which has been adopted is 
in agreement with the above definition of the telluroid. 


Dermanis (1984) has rigorously derived the fundamental boundary 
condition on the telluroid for the Marussi mapping in terms of 
curvature parameters [Dermanis, 1984, eq. 53]. He has also given the 
corresponding form of the boundary condition on the reference 
ellipsoid in terms of ellipsoidal, geodetic and spherical coordinates 
[ibid, eq. 82, 85, 86]. 


2.2 Gravity Anomaly and Potential Coefficient Relations 

The linearization procedure as presented before has lead to the 
boundary condition 


AW + c 


p 


with 




(2.44) 


Ag = Ifpl - I?qI 
AW = W 0 - U 0 

r,, aw . aw aw 1 

e P - [(l cos e) 3h f M3+ v Ncos+axJp 
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(2.45) 


valid under the approximation (2.24) on the telluroid (E) as defined by 
(2.29), (2.40), and (2.41). The geometry of the original problem (Figure 
1) has now been simplified, due to the approximations (2.24) and (2.29), 
and is now described in Figure 3. 


However, since both the geometrical shape and the gravity field of 
the reference ellipsoid differ from the corresponding features of a 
homogeneous rotating ball by quantities on the order of the flattening 
(~ 1/300), one can expand geometrical and dynamcial quantities refering 
to the ellipsoid around their "spherical" values in a Taylor series in 
terms of any parameter characterizing the departures from the 
"spherical case" (e.g. f, e a , e' 2 , etc.). Truncation of such series to 



Figure 3. The Geometry Associated with the Boundary Condition (2.44). 


21 

different orders, results in different approximations. For example, 
omission of all but the zeroth term results in the well known "spherical 
approximation" [Moritz, 1980, Sec. 4.2]. It should be carefully kept in 
mind that such series expansions are not meant to change either the 
reference surface, which always remains the telluroid, or the 
computation of normal gravity at the telluroid points. Such 

computations should be carried out to a high degree of accuracy. 

The importance of these series expansions rests on the fact that 
their zeroth order (spherical) terms provide, through the boundary 
condition, a simple way of relating gravity anomalies to potential 
coefficients through a spherical harmonic representation. The 
remaining ellipsoidal terms due to their small magnitude can be treated 
otherwise. 

In the following, the form of the boundary condition when 
ellipsoidal terms to 0(e a ) are modeled is presented, following Jekeli’s 
(1981) formulation. 


2.2.1 The Boundary Condition Considering Ellipsoidal Terms to Order 
e* 

Consider through the telluroid point Q a coordinate surface 
[Heiskanen and Moritz, 1967, Sec. 1-19], u = b 0 (Figure 4). For each 
telluroid point Q<, the corresponding coordinate surface will have 
constant linear eccentricity (by definition), but a different semi-minor 
axis uj = b Q1 . Accordingly, since 

E = b 0l (l-ej)"*e 1 (2.46) 

each coordinate surface will be characterized by a different 
eccentricity e^ However, to the accuracies involved here and for 
points close to the earth’s surface, the eccentricity of all coordinate 
surfaces can be considered constant, equal to the eccentricity of the 
reference ellipsoid [Jekeli, 1981, p. 122]. Consider a local spherical 
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coordinate system at Q (Figure 4) with unit vectors <t r{j in the 
direction of the geocentric radius and positive outwards, ^Q^along the 



Figure 4. The Coordinate Surface Through the Telluroid Point Q. 


increasing geocentric co-latitude and along the increasing 


longitude. Then 

^h Q = cosV- Q ^r Q - sin (2.47) 
and the gradient vector with respect to this system becomes 

( 8fad •)« = (&),*.-, + (^i) *8 a + (2 ' 48) 

Since, 

(fs), = (Srad •)„ ■ i h(> (2.49) 

we have 

(Ik),, = cosl, « (l?] a - aiTr *° (?f§) a (2 - 50 > 


It can also be shown easily that to the order of e 2 , 



sin^Q 8 e 2 sin0QCos0 Q 
cosY'q * 1 
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(2.51) 


Accordingly, 

(!k) q = (fe) Q ‘ ea sin0 Q cos0 Q [^} q + 0(e 4 ) (2.52) 

On the other hand, the normal gravity potential U is given by 
[Heiskanen and Moritz, 1967, Sec. 2-9] 

U(r, 0) = ^-[l - J^Janf^'p^tcos©)] + ± w 2 r 2 sin 2 0 (2.53) 

where G is the Universal gravitational constant 

M' is the mass of the reference ellipsoid (r.e.) 
a is the semi-major axis of the r.e. 
u is the rotational rate of the r.e. 

J 2n are the even degree zonal harmonic coefficients of the 
normal gravity field 

P 2n (cos0) are the even degree Legendre polynomials. 


However, 

becomes 

since J 2n - 0(e 2r> ), equation (2.53) truncated to 

the 0(e 2 ) 

U(r, 

0 ) ~ 

^-[l - J*(~j P 2 (cos0)j + i W 2 r 2 sin 2 0 

(2.54) 

Hence, 




7 - ~ 

au 

‘ ah ' 

= - + e 2 sin0cos0 yields for point Q 


„ GM' 

7Q ~ ^7 

[i- 

3J 2 [— -] P 2 (cos0(j)j - w 2 rQsin 2 0Q (1 - e 2 cos 2 0 Q ) 

(2.55) 


and 

- e 2 sin0cos0 yields for Q 


(fh) ~ ~ ~ 6J 2 (^“) P 2 (cos 0 q )] - <v 2 sin 2 0 o (l - 3e 2 cos 2 0 Q ) (2.56) 



Inserting (2.55) and (2.56) to (2.44) and carrying out the algebra 
consistently omitting terms smaller than 0(e 2 ), one finally gets 



+ e 2 sine Q c OS e Q (|y q + [6J 2 ^ P 2 (cos0 Q ) - sin 2 0 Q ]T Q 

+ ~ AW - [6J a ^ P a (cos0 Q ) - sin 2 0 Q ] AW + c P (2.57) 

which is the form of the fundamental boundary condition when 
ellipsoidal terms are considered only up to 0(e 2 ). The following 
abbreviations will be used: 


| *h = e 2 sin©cos0 (^||) (a) 

ey = [6J 2 ~ P 2 (cos0) - sin 2 ©] T (b) 

i 

i 


(2.58) 


2.2.2 Spherical Harmonic Representation of the Gravity Anomaly 

For the moment let us assume that there are no masses outside of 
the telluroid. Then, if V and V' are the gravitational potentials of the 
earth and the reference ellipsoid respectively, then they are both 
harmonic in the space outside of the telluroid. Hence they can both be 
expanded into series of solid spherical harmonics [Heiskanen and 
Moritz, 1967, Sec. 1-9] as 


V(r,8,X) . ^ * jjf 

;] X (C nm cosmX + S nm sinrnA)P nm (cos0)] 

m=0 J 

(2.59) 

and 



v< r . 0 > = 2f[i + l { 

r 1 n=l l 

a] 2 "-. - 1 

7J C 2n P 2n (cos0)J 

(2.60) 


where 
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is the mass of the earth (in contrast to M' which denotes the 
mass of the ref. ellipsoid) 

are the fully normalized associate Legendre functions of the 
first kind 

are the fully normalized potential coefficients (referring to 
the earth's gravitational field) 

are the corresponding coefficients of the ellipsoidal field 
always denote geocentric radius, co-latitude and longitude 
respectively. 

The difference 

AM = M - M' (2.61) 

is assumed to be sufficiently small that for the first and higher order 
terms M' can be replaced by M in (2.60). 

If the rotational rate of the reference ellipsoid is assumed to be 
the same as the actual rotational rate of the earth then 

T=W-U=V-V' (2.62) 

As the difference of two harmonic functions is harmonic in the space 
outside (2) so that 

T(r,0,X) = ~ Z (f] n I (C* w cosmX + S nni sinmX)P nni (cos0) (2.63) 

r r n=t m=o 

with 

f C nro - c' m if m = 0, n = 2k, k e N 
CL = I (2.64) 

l C nm otherwise 

Inserting (2.63) into (2.57) one gets 

Ag = - fnr + ~a 2 (n-i) (~~) ° I (C^cosmAq + S nro sinmX Q )P nm (cos0 Q ) 

r Q r Q n=2 kt Q j m-o 

+ ( c h) Q + (* 7 )« + f “ AW-[6J 2 ^ P 2 (cos0q) - 3 “^ q2 sin 2 0q] AW 
p 


M 



(r,0,X) 


+ e 


(2.65) 
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where M' has been replaced by M in the sixth term following similar 
reasoning as before (equation (2.60)). 

2.3 The Mathematical Model - Practical Aspects 

Equation (2.65) constitutes the fundamental mathematical model 
which establishes the functional relationship between the gravity 
anomaly Ag and the potential coefficients (C* m , S nm ). Therefore, it 
provides the basis for the formulation of observation equations to be 
used in a least squares adjustment for the estimation of {C* m , S nm ) 
from A g. The variables appearing in (2.65) are classified as follows: 

(a) Observed quantities are considered the gravity anomalies Ag, while 
unknown parameters are the potential coefficients {C* m , S nm ), as well as 
the quantities GAM = GM - GM' and AW = W c - U 0 . 

(b) The quantities {a, GM, J 2 , w) are constants and refer to the 
adopted mean earth ellipsoid, while {r, Q, \}q are considered known for 
each telluroid point Q. 

Prior to the formation of the observation equation it is necessary 
to examine how the available data comply with the requirements and 
assumptions of the mathematical model and to determine the reductions 
that need to be applied to the data and the modifications that should 
be made to the model so that any incompatibilities between the model 
and the data are removed. 

Also, the varying magnitude and information content of the 
components of the mathematical model justify special treatment of 
certain terms (such as c h and e^), so that a more computationally 
convenient formulation can be achieved. 

These aspects will be discussed in detail in the following sections. 
As a result of this discussion the final formulation that was used in 
the computations will be presented. 
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2.3.1 The Observable 

The primary observable from terrestrial measurements is the 
magnitude of the gravity acceleration g P = l^pl at the surface point P 
(see Figure 5). 

To evaluate the surface free— air gravity anomaly Ag [Heiskanen and 
Moritz, 1967, p. 293], which appears in equation (2.65), it is needed to 
evaluate the magnitude of the normal gravity acceleration y Q = 1% 1 at 
the corresponding teliuroid point Q. 

If 7 q o denotes the normal gravity at the footpoint Q 0 on the 
surface of the reference ellipsoid (Figure 5), then 



Figure 5. Geometry Associated With the Observable. 
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( 2 . 66 ) 


since the distance IQ 0 dl is, according to the definition used here 
for the telluroid, the normal height H* [ibid, p. 292] of the surface 
point P. However, it is known [ibid, p. 328], that H* (which is usually 
unknown) differs from the orthometric height H at P (which is more 
readily available) by about one meter in the worst cases (mountainous 
areas). Accordingly, without significant loss of accuracy one has 


7q = 7q q + 



(2.67) 


and the surface free-air anomaly becomes 


Ag = g p 



or [ibid, p. 293, eq. 8-9] 


( 2 . 68 ) 


Ag = g p - 7Q 0 [l - 2(1 + f + m - 2fsin 2 * Qo ) |] - 3-y Q(> (g) 2 (2.69) 

where 

f = (2.70) 


2 a 2 b 

GM' 


and 7 q q 


can be evaluated using Somigliana’s closed formula 
a7 e cos 2 +Q + bypSin 2 <>Q 
(a 2 cos 2 +Q + b 2 sin 2 +Q) 1 / 2 


(2.71) 


(2.72) 


with the equatorial and polar normal gravity y e and y p given by 

r.*S (l---!*?) <2.73, 


_ a f, + a sUil 

r P a 2 t 3 q 0 J 


(2.74) 
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with [Moritz, 1984, pp. 388-398] 

q ° = 1 [(* + ^sr) arctan e * “ I 7 "] (2.75) 

qo = 3(l + 4r)(l - ~r arctan e'J - 1. (2.76) 

Apart from the approximation of H* by H, evaluation of the surface 
free-air gravity anomaly by equation (2.69) would be in agreement with 
the requirement of the mathematical model (2.65). 

Additional discussion related to the surface free-air anomaly will be 
made in Chapter III where the description of the actual data used in 
this study is given. 


2.3.2 Atmospheric Effect 

In deriving equation (2.65) it was assumed that there are no 
masses outside the telluroid. However, even if we neglect the fact that 
the telluroid can be below the surface of the earth, there still exists 
the effect of the atmosphere that needs to be taken into account. 

To illustrate the effect we will use a simplified atmospheric model 
which assumes spherical stratification [Moritz, 1980, pp. 422-425]. 



Figure 6. Atmospheric Correction. 
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The magnitude of the atmospheric attraction at the surface point P is 
in general (Figure 6) 

Fat* = G (2-77) 

F P 

where m(r P ) is the atmospheric mass below the sphere passing through 
P (radius r p ). 

If M a denotes the total mass of the atmosphere then 

Fat* = G 7*1 - G (2.78) 

r P r P 

where M(r p ) is the atmospheric mass above the sphere of radius r p . 

F a t m needs to be subtracted from the observed gravity g p for the 
effect of the atmosphere to be removed. However, if in the definition 
of the normal potential, the mass of the reference ellipsoid includes the 
mass of the atmosphere (as in GRS67 and GRS80) then, to remove the 
atmospheric effect from the gravity anomaly Ag it is required that the 
term 

= G (2.79) 

r P 

be added to the gravity anomaly as obtained from (2.69). 

Tabular values of Sg A are given, as a function of the elevation 
above sea level [IAG, 1971, p. 72]. These values have been computed 
based on a nearly-ellipsoidal stratification model and the use of two 
standard atmospheres, the CIRA 1961 and the U.S. Standard 
Atmosphere. The tabulated values represent rounded averages from 
these two atmospheric models. For computer implementation a 
convenient way to interpolate is to use a polynomial of the elevation 
that fits the tabulated values. Wichiencharoen (1982) performed such a 
polynomial fit; the following quadratic is the result; 

<5«A(mgais) = 0.8658 - 9 . 727xlO' 5 H( m ) + 3.482xl0“ 9 H 2 ( m ) . (2.80) 

This quadratic has been used in this study for the calculation of the 
atmospheric corrections 5g A * As it is obvious from (2.80) <5g A amounts 
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to about 0.9 mgals at sea level (where the correction is maximum). The 
indirect effect (shifting of the level surfaces due to the condensation 
of the atmospheric mass on the earth’s surface) will not be considered 
due to its small magnitude (-0.7 cm at sea level) [Moritz, 1980 p. 425]. 

2.3.3 Ellipsoidal Effects 

As it can be seen from equations (2.58a, b) the unknown potential 
coefficients (Cj[ m , S nm ) appear on the right side of equation (2.65) not 
only in the harmonic sum but also in the terms e h and tj. However, 
the magnitude and information content (see also Appendix A) of these 
terms suggest that it is neither reasonable from the estimability point 
of view, nor computationally convenient, for such terms to be 
incorporated into the design matrix of the adjustment. Rather, it is 
preferable to evaluate and ty beforehand, based on some prior 
information for the unknown coefficients and apply these terms as 
corrections to the observations. 

It will be assumed here that a set of potential coefficients complete 
to maximum degree and order N max is known (e.g. Rapp81 or OSU86F). 
Assuming, for the present purpose that these coefficients refer to a 
mean earth ellipsoid with mass equal to the mass of the earth and 
center coinciding with the center of mass of the earth one has (see 
also equation (2.63)). 

T(r,8,X) = — Z - Z (C* m cosmX + S nm sinmX)P nm (cos0) (2.81) 

r n=2 '•r.* m=o 

where (~ ) is used to denote that the truncated series with the a priori 
coefficients provides only an approximation to the true value of the 
disturbing potential. Using the well known relations 

Pao(cos8) = | cos a 0 - h (a) 

(2.82) 

sin a 0 = 1 - cos a 0 (b) 


and the abbreviations 
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e 2 GM 

c i - r 2 

(a) 

c 2 = 3J 2 r4 + u 2 r 

(b) 

c 3 = Ja + " 2 r 

(c) 


(2.83) 


one has from (2.58a) and (2.58b) 




Nmaxf a ) n n - 

= c x E - Z (C£ m cosmA + S nm sinmA)sin0cos0 

n=2 lr; m— o 


dP nm (cos0) 


d0 


(2.84) 


and 


[ Nmaxfo'in n - - 

C 2 l - l (C* m cosmA + S nm sinmA)cos 2 0P nm (cos0) 

n = 2 ^ r; m-o 

Nmaxfoln n - . - - 1 , _ s 

-c 3 l - I (C* m cosmA + S nm sinmA)P nm (cos0) j (2.85) 

n = 2^ m=o J 

Using the normalization factors 

„ _ / (2n+l)(n-m-l)(n-mT . „ ^ 

nm (2n-3) (n+m-1) (n+m) * 

r ( 2n+l ) ( n+m+1 ) ( n+m+2 ) . n . 

v nm " J (2n+5) (n-m+1) (n-m+2) n , m n 


( 2 . 86 ) 


dP 


(2.87) 

and the recursive relations for sin0cos0 [Moritz, 1980, p. 321, 

eq. 39-46] and for cos 2 0P nm [ibid, p. 326, eq. 39-76] one finally has 

dP„ m _ 


sin0COS0 - 8 nm V nm P n + 2 f m ^ c nm’*nnJ > n - i,i 

d0 


( 2 . 88 ) 


with [ibid, p. 321, eq. 39-47] 


n ( n-m+1 ) ( n-m+2 ) 
(2n+l) (2n+3) 


n 2 -3m 2 +n 

nm = (2n+3) (2n-l) 

(n+1) (n+m) (n+m-1) 
"« _ (2n+l) (2n-l) 


(2.89) 
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and 

COS 2 0P nm - <X nn ,V nm Pn+a , m + ^nm^nra + 7’nm u nmPn-2 , m 

with [ibid, p. 326, eq. 39-77] 

(n-m+1) (n-m+2) 

(2n+l)(2n+3) 

n _ 2n 2 -2m 2 +2n-l 
Pnm (2n+3)(2n-l) 

(n+m) (n+m-1) 

rnm (2n+l) (2n-l) J. 

Using the abbreviation 

®nm(^) = C* ro cosmX + S nm sinmX 
from (2.84) and (2.88) one gets 

[ Nmaxfa'jn n * 

£ irJ ^ ®nm(^) a nni v nmPn+2 J in 

n=2 KrJ m=o * 


Nmax n n * - 

+ Z f Z Dnn,(X)b nm P nm 

n=2 vrJ m=o 


Nmaxfoln n~ 2* _ i 

+ Z If J Z D nm (X)c nm u nm P n _ 2 J 

m -2 m=o ’ J 

while from (2.85) and (2.90) one gets 

( rNmaxfp)n n ~ 

c 2 [ X ["J £ D nm ( X)c( nm V nn ,Pn+2 , m 

l n =2 KrJ m- o * 

Nmaxfoln n - 

+ Z If I D nm (X)^ nm P nm 

n - 2 vrJ m=o 

Nma x f al n n“2 ~ - 1 

+ ^ irJ ^ ®nm(^)7nm u nmPn-2 , m I 

n=2 m=0 * J 

Nmax fain n ~ - 1 

- c 3 l f l D nm (X)P n J 

Using equations (2.93) and (2.94), the ellipsoidal terms s h 
be computed from a set of potential coefficients complete 


(2.90) 


(2.91) 


(2.92) 


(2.93) 


(2.94) 

and can 
to maximum 
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degree and order N max . 

As it will be seen in the next section, the actual evaluation of the 
ellipsoidal terms e h and ty in this study was made in terms of area 
mean values (see equations (2.111) and (2.112)), so that these 
corrections are compatible with the data used here. 

From the previous discussion it becomes clear that the approach 
adopted here for the computation of the ellipsoidal corrections is the 
reverse of the one proposed by Cruz (1986) and implemented in the 
development of the OSU86 series of high degree geopotential models 
[Rapp and Cruz, 1986a, b]. In this formulation the corrections are 
applied to the gravity data and not to the potential coefficients 
obtained from them. The reason, being rather obvious, is that the 
combination solutions for the TOPEX gravity models are to be 
performed at NASA/ (Goddard Space Flight Center) and not at the Ohio 
State University. Hence, operationally corrections to the coefficients 
obtained from these solutions cannot be applied. The normal equations 
formed from surface gravity should include any systematic reductions 
necessary to make them compatible with the satellite data derived 
normals. 

In Appendix A numerical values related to the magnitude of these 
corrections and Figures illustrating their geographical distribution are 
given. These, compared to the corresponding Figures given in [Cruz, 
1986] verify the equivalence between the two methods. 


2.3.4 Transition from Point to Area Mean Values 

To this point the discussion on the modeling of surface gravity 
data was restricted to point value considerations. The gravity anomaly 
4g in equation (2.65) represents a point value referring to the telluroid 
point Q or equivalently (see discussion after equation (2.13)) to a 
surface point P. However, for global gravity modeling one is 
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restricted, in practice, to use area-mean values of Ag, due to the lack 
of point data and economy considerations. 

Area-mean values of Ag, hereon to be denoted by Ag, are usually 
defined over either equi-angular or equal-area blocks. In view of the 
data to be used in this study (see also Chapter III), it will be 
considered in the following that Ag refers to equi-angular blocks, i.e. 
blocks bounded by meridional and parallel arcs of equal angular extent. 

Consider an equi-angular grid on the ellipsoid and let AX be the 
angular distance between two successive meridians. Such a grid forms 
N = tt/AX latitudinal bands and 2N longitudinal sectors on the ellipsoid, 
creating in total 2N 2 equi-angular blocks. Each block will be hereon 
identified by two subscripts i and j, implicitly defining the latitude 
and longitude bounds of the block, so that ♦o - w/2 , = -tt/ 2 and X 0 

— 0, X 2 ^ — 2rr. 

Let H { j denote the mean orthometric height of the ij th block. 
Then Ag<j is interpreted to represent the average behavior (value) of 
the spatial function Ag(r ,0,X) over the surface <s \ j located at constant 
distance H ( j from the surface of the reference ellipsoid (Figure 7). 



‘ j* 


Figure 7 Geometry of the Surface Element <s 
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It is obvious that although Hjj is considered constant over the 
block, the geocentric distance r varies (in the latitudinal direction) 
over the surface element <Xij. Let 


♦ i = ^ (♦i + +i+i ) (a) 
Xj = 2 (^j + ^j+i) (b) 


(2.95) 


be the geodetic coordinates of the center of the ij th block. Then 


r <j 

where 


<*ij 




z u 


1/2 


(2.96) 


Xjj = (N i + j)cos<^jCOsXj 

Yj j = (N< + Hi j )cos+iSinXj 

Zi j = [Ni(l-e 2 ) + Hijlsin^i 

For an arbitrary point on with geodetic latitude ♦ we have 


(2.97) 


r(4») = ?ij 



, 


(♦ - ♦i ) + . . . 


(2.98) 


One can easily find [see also Rapp, 1983b, p. 67, eq. 48], that 

*■(♦) = r u - | e2s ^ i Ni (M, + H< j ) (♦ - ♦,) + ... (2.99) 

with 

_ a 

N, = (2.100) 

(l-e 2 sin 2 *i ) »/ 2 

a(l-e 2 ) 

Mi = (2.101) 

(l-e 2 sin 2 4»i) 3 / 2 

The term 3r/3$ as it can be seen from (2.99) becomes maximum (in the 
absolute sense) around ♦=45”, where its value is about -370 m/deg. 
Here the variation of r over the surface <X{j will be disregarded and it 
will be assumed that r = r ( j = const, over <r s j . Additional discussion 
on this approximation is given in Appendix A. 
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Now consider equation (2.65) and apply on both sides the integral 
operator 

Int {•} = JJ* sin0d0dX (2.102) 


which represents the surface integral over the surface element a { j 
(where r = Fj j = const.). If (2.65) is divided by the spherical area of 
<r | j (given by Int{l}), then the left side becomes 


Int{Ag(r=r, j, 0, X)} 
Int{l) 


(2.103) 


which represents the area-mean value Ag { j of Ag(r ,0,X) over the 
surface element <r,j. Accordingly, from (2.65) one gets 


= h~ r~* 2 (n-i)(^-) n i (c* m ic i + s n „isj)ipj, 

J r ij n = 2 m=o 

+ IE^ J + IE.J, j + aiGfiM + a 2 AW + IEp j 


where 


rr (&i+i f A j +l 

Affj = Int(l) = JJsin0d0dX = J sin0d0 j dX 

a \ Xj 



mXdX 


IPnm = P nm (cos0)sin0d0 


< J = 4^7 Int W> 



(2.104) 

(2.105) 

(2.106) 

(2.107) 

(2.108) 


(2.109) 
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a 


2 





(cos 2 0j 


+ COsOjCOsSj+x + cos 2 0 j + i) 


+ 3 J 2 



cj2r i j a 
GM , 


(2.110) 


The terms IE^J and IE^J are computed from equations (2*93) and (2.94) 
respectively, after averaging over the surface element j* However, it 
can be verified numerically that without significant loss of accuracy 
these terms may be evaluated on the surface of the ellipsoid rather 
than on the ground level* Hence 


j 1 / rNmaxf « n n * -4 

IE; J = c t I hr- Z ID„ m a nm v nm IP „+a m 

A<r i L n = 2 m=0 9 

Nmax ( o 1 n n - \ - 4 

+ l hr- Z iDJ m b nm ip; m 

n=2 ^ r E i ; m=o 

Nmax f n 'in n ~ 2 - 4 1 

+ l hr- l iDLc nm u nm iP,;- 2 , m j (2.111) 

n — 2 1 E j m— o 

*44 1 f / fNmaxf n In n *4 -4 

IE r - IT" 3 c 4 ^ Z IDnm^nmVnmIPn+2, m 

' 1 l n = 2 ^ r E j ‘ <n-o ’ 

Nmax ( q 1 n n ~ 4 - 4 

+ Z p 2 - Z IDnm^nmlPnro 

n = 2 Vr E , ; m=o 

Nmax f n 'In n — 2 * 4 - 4 *| 

+ I p 2 - l IDLrnn,U„ m IPn-2, ffl 

n = 2 '• r E, ; m-o J 


- Nmax falnn ~ 4 -il 

- c 3 Z M- Z I^nmlPnin j 

n = 2 ' r E j m=0 1 


where 


ID 


dL = c* icA + 


c’ ICi + S n _IS„ 

^ n m m n m -^^m 


( 2 . 112 ) 


(2.113) 


and r Ei is obtained from (2.96) and (2.97) for H { j =0. The terms c(, 
c 2 , c 3 are computed from (2.83) by replacing r with r E i* The 
geocentric co-latitude 9 appearing in the above equations is computed 
from the geodetic latitude by 
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d = ^ ~ arctan[(l-e a )tan+] (2.114) 

In addition, since we have assumed constant orthometric height Hj j 
over the ij th block, the atmospheric correction 6 gj( J to be applied to 
Ag ( j will be given by equation (2.80) for H = Hjj. Hence, the 
mathematical model relating area-mean values of surface free-air 
anomalies to potential coefficients becomes 


Ag {j + 5gi J = ^7 A I (n-1) (**-)" I (C*JC 1 + S nm IS^)IPn, 

fl<r i r 1J n = 2 ^ r 1j' m=o 


+ IEi J + IE^ j + aiGtfM + a 2 AW + IEp J 


(2.115) 


2.3.5 Vertical Datum Inconsistencies 

Modeling of a global set of gravity anomalies based on equation 
(2.115) assumes that all gravity anomalies Ag ( j refer to a unique global 
vertical datum. In the absence of such a datum it is possible that 
subsets of a global gravity data base contain anomalies referring to 
different (regional) vertical datums. Such inconsistencies between the 
data have to be considered in the mathematical modeling [Rapp, 1984]. 

Let H* be the normal height [Heiskanen and Moritz, 1967, p. 292] of 
a point P with respect to a unique (but unknown) fundamental level 
surface, and the normal height of the same point referring to the 
k th regional level surface. Then the normal gravity jq for the 
corresponding telluroid point Q is (see equation ( 2 . 66 )) 
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+ [(IS), + 5 (0), 

M o 

or 

r. = r& + [&) a - (0), (2.116) 

“o o 

where y|j is the normal gravity computed from elevation information 
referring to the regional level surface* Hence, 

a*’ 1 *" - [(??), + (IS) a (2.117) 

M o M o 

where Ag k denotes the gravity anomaly that refers to the k th regional 
reference level surface. Parametrization in terms of potential 

differences is also possible since [Heiskanen and Moritz, 1967, eq. 4-44] 


Cp 

H* = 1+ (1+f+m- 2fsin a $) 

L 

CP + ( CP 11 
a 7Q 0 + a 7 Q o J 

(2.118) 

and 



C p = W 0 - W p 


(2.119) 

/"i k _ tj k _ rj 

Cp - »o w p 


(2.120) 

so that 



ACj = C p - Cp = W 0 - Wo = AWo 


(2.121) 

Accordingly, equation (2.117) becomes 



ig » Ag k - [(J$ a + (0) a Hj) 

AWo 

(2.122) 


In terms of area-mean values equation (2.122) takes the form 
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Agjj - Agjj + 


~ 3J 2 f"3 + 


2 « a ?E? 

a 2 1 


r E % GM 


(cos 2 0 1 +cos 0 i cos 0 i + l +cos 2 0 1 + 1 ) 


a 2 " 2?e i 

+ + -®r 


- 4 \ 


AWo 


(2.123) 


According to (2.123) equation (2.115) becomes 


Agfj + figl J = I (n-1) («*-)" I (C* m IcJ + S nro IsJ)IP;, 

J “i r lj n = 2 lr 1j' «n=o 


+ IEh J + IE~ J + axGfiM + a 2 AW + a 3 Aw£ + Ie1 J (2.124) 


with 


a 3 = - ■= — + |3J 2 ^—3 + 


,2 " 2? E ai 


GM 


(COS 2 ©! + COSQ^OS©^! + COS 2 ©^) 




- 3|J 2 


cj‘r E 

“gm" 


+ 6 S| 


(2.125) 


As it can be seen from equations (2.110) and (2.125) 


a 2 ** -a 3 


(2.126) 


a relationship which expresses none other than the fact that AW and 
AWfc cannot be recovered simultaneously from gravity anomaly 
observations alone. Rather their linear combination 


AW - AWo = Wo - U 0 


(2.127) 


is the only parameter that one might attempt to recover from such data. 
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2.3.6 Computational Formulas 

Equation (2.124) represents the mathematical relationship between 
surface mean free-air gravity anomalies and potential coefficients, when 
atmospheric effects, ellipsoidal effects to 0(e 2 ) and vertical datum 
inconsistencies are taken into account. It is valid with respect to an 
arbitrary reference ellipsoid since the terms G<5M and AW are accounted 
for. Before presenting the final formulation used in this study two 
specific terms in (2.124) must be considered. 

The term + a 2 AW will give rise to a nearly constant term in 

the gravity anomaly if evaluation of (2.124) is made on the surface of a 
sphere (radius R). It can be easily seen that in such case 

a x GAM + a 2 AW = |~ + | AW -[6J 2 IP 2Q - IP 22 ]aW. (2.128) 

The values of GAM and AW are anticipated to be in the order of the 
uncertainties of the currently used best estimates of these parameters 
(see also Rapp, 1983a). Consequently, the very long wavelength 
latitudinal variation caused by the last term can be neglected in view 
also of the magnitude of its coefficient. Note that the average gravity 
anomaly over the sphere R is independent of this term, since its effect 
will be eliminated if a x GAM + a 2 AW is averaged over the whole sphere. 

If the latitude dependence of a x GAM + a 2 AW is neglected, one has 

a x GAM + a 2 AW * - + JL AW 

r ’ j r 'i 


or 


a x GAM + a 2 AW * 


GM f GAM ] 
? * j 2 l GM J 


GM [ 2RAW ] 
r,jH l GM J 


(2.129) 


In the analysis of gravity anomalies located on the surface of the 
earth, the difference between Fjj and R is disregarded and (2.129) is 
written as 


a,G<5M + a 2 AW « - 
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GAM _ 2RAW ) 
GM GM J 


(2.130) 


with R being a mean earth radius. This term can now be incorporated 
in the harmonic summation in equation (2.124) so that 


h k ii + 


1 GM 
A<r, r,j 2 



+ s nm isJ,)ip; 


+ l£j J + I&J, j + a 3 AWo + IEj j (2.131) 

where 


-* _ GAM _ 2RAW 
00 " GM GM ’ 


(2.132) 


If AW=0, then equation (2.131) is rigorously equivalent to (2.124). 
However, if AW*0, then (2.131) and (2.132) provide only an 

approximation to (2.124). As it is explained later, in the analysis of an 
incomplete set of gravity anomalies on the surface of the earth it is 
crucial to incorporate the C? 0 term in the modeling, even if G<5M=AW=0, 
due to the correlations between the unknowns arising from the 
incomplete coverage. The recovered C* 0 from such analysis may then 
be interpreted according to (2.132), mainly to check the consistency of 
the solution, since physical constants such as GM and W 0 (which is 
related to the semi-major axis a) are more accurately determined from 
satellite techniques [Rapp, 1983a]. 


On the other hand, the term * p is primarily a function of the 
deflection of the vertical as it can be seen from equation (2.34). Prom 
[Jekeli, 1981, eq. 4.22, with the sign of the fourth term corrected] one 
has 


,„=■ (2.133) 

since W=U+T and U is longitude independent. To examine the 
magnitude and information content of e p , it is helpful to consider the 
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following approximations 


St * -g 

ah 8 


(a) 

3T 0 
M3* 

aT 

rae 

(b) 

3T 

NCOS* ax ~ 7V 


(c) 

au au .. 

M3* _ M30 ~ 

au 

rae 

(d) 

1 - cos 8 ~ — 

e 2 

(e) 


According to the above, equation (2.133) becomes 


The first term is (2,135) is approximately equal to 
0 2 =£ 2 +tj 2 , and g p s 7 p). With gp s 7p*9.8 m/sec 2 and 
magnitude of this term is found to be about 2 pgals, so 
term can be safely neglected. 


On the other hand, from equation (2.54) one has 
|f * Ja (“) + u 2 r 2 ]sin0cos0 

and using (2.55) one gets 

1 3D„ r_ T a 2 , cj a r 4 l . . _ 

— -rrr ~ 3J 2 — + sin0cos0. 

y 36 L r GM J 

Accordingly, equation (2.135) becomes 

£ p * {[ 3J2 (?) 2 + l ^] sin0cos0 Help . 


(2.134) 


(2.135) 

9 2 g/2 (since 
£=77=10", the 
that such a 


(2.136) 


(2.137) 


(2.138) 
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Using a nominal value of 6371 km for r one can verify numerically that 


3J 


.(?)* 


cj 2 r 3 

GM' 


(2.139) 


Comparing now (2.58a) to (2.138) (in view also of (2.139)), one concludes 
that the effects of e h and e p are almost identical both in terms of 
frequency content and magnitude. Area-mean values of c p can thus be 
computed by rescaling the area-mean values of e h , i.e. 


If 



where 


(2.140) 


Cl 


3J 2 


a*GM 


CJ 


2 ?e 




and c[ is given by (2.83a) with r = r Ej . 


(2.141) 


Summarizing, the mathematical model relating area-mean values of 
surface free-air anomalies to potential coefficients, considered in this 
study is 


AtflJ 


_ 1 GM 2 , a 1" a T „J 

g,J " Aff | r f j 2 „| ? (n 1 If j j) Jo (Cni " IC -" 


J nm J 


)IPn, 


+ a 3 AWo 


(2.142) 


where <5g?j represents the systematic reductions necessary to be 
applied to the surface mean free-air anomalies and is given by 


<5g?j = <5gi j - (IE h + IE 7 + IE p ) 1J or 

««?, = Hi 1 - [(l + ff)l8„ + re,]] 1 


(2.143) 
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A step-by-step procedure for the numerical evaluation of the terms in 
(2.142) and (2.143) can be outlined as follows; (the assumption is made 
in the following that an approximate set of coefficients complete to 
degree and order N max , and a set of area-mean values of orthometric 
heights are available). 


(i) Evaluate IC^ and IS„ by 


ICi 


2 . mAX r 

— sin cosmX 

m 2 

AX 


J 


and 


if m * 0 
if m = 0 


(2.144) 


IS l = 

(ii) Evaluate IPj, m . This can be done efficiently using the recurrence 
relations of Paul (1978). Note that IPj 0 yields the spherical area 
A<Xj of the ij th block. 

(iii) Evaluate r S j from equations (2.95) through (2.97). From the same 
equations, with H f j =0, evaluate r E < . 


2 . mAX . 7 

- sm — simA 
m z 

0 


J 


if in * 0 
if m = 0 


(2.145) 


(iv) Evaluate <5gj(j from (2.80) and IE^J, IE^J from equations (2.111) 
and (2.112). Consequently, the systematic reduction <5gfj to the 
mean anomaly can be evaluated by (2.143). 

(v) If vertical datum inconsistencies are to be parameterized then a 3 
in (2.142) can be evaluated from equation (2.125). 


A setup as the one described above will enable one to numerically 
evaluate all coefficients necessary for the formation of observation 
equations based on the model (2.142). 
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Concluding this section an important remark should be made about 
equation (2.142). Theoretically, the convergence of the infinite series 
in (2.142) is guaranteed only in the region outside the smallest sphere 
enclosing all the generating masses of the potential (Brillouin sphere). 
However, to now no successful theoretical argument has been given 
proving or disproving the convergence of the series on the earth’s 
surface [Jekeli, 1981, p.4]. The intention in this study is to use a 
truncated version of (2.142) (to N max a 50), for evaluation on the 
earth’s surface. To this extent the problem of convergence will not be 
addressed here in view also of the numerical results reported by Jekeli 
(1981). 


2.4 The Alternative Use of Terrain Corrected Anomalies 

It was proposed by Moritz (1966, p. 104) that linear solutions to 
the bvp of Molodensky can be obtained using anomalies reduced to the 
surface of the reference ellipsoid through the free-air gradient, as 
follows 

Ag* = Ag - h (2.146) 

where Ag is the surface free-air anomaly and Ag* the downward 
continued anomaly. 

Moritz [ibid, p. 105, equation 280] has proved that 

Ag* = Ag + Gj + | l n(hAg)„ (2.147) 

“ n=o 

where Gi is the Molodensky correction term [Heiskanen and Moritz, 
1967, p. 307]. 


For lower degree harmonics one can neglect the term 
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s l n(hAg) n (2.148) 

K n = 0 

as being of the order of h/R. Hence, from (2.147) one obtains 

Ag* = Ag + G x (2.149) 

On the other hand, it is shown by Moritz (1966, p. 105, equation 285b) 
that 


G x = G' - l n(AhAg) n (2.150) 

6tl n = 0 

where 

Ah = h - h m (2.151) 

and h m is a mean elevation of the area in question. G' is a correction 
term associated with the "Pellinen type solutions" to the bvp of 
Molodensky [ibid, p. 92], Neglecting the term 

- | o n(AhAg) n (2.152) 

with the same reasoning as in (2.148), for lower degree harmonics, one 
has 


Ag* = Ag + G' (2.153) 

or in view of (2.146) 

-|££h = G\ (2.154) 


In addition, if the assumption is made that the surface free-air 
anomalies are linearly correlated with elevation, i.e. 
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Ag — a + bh (2.155) 

with a,b being approximately constant, then it is proven in [ibid, 
equation 262] that G' becomes identical to the conventional terrain 
correction tc [Heiskanen and Moritz, 1967, p. 131]. Hence for lower 
degree harmonics and under the assumption (2.155) one has 

Ag* — Ag + tc. (2.156) 

Accordingly, if terrain corrections tcjj are available in terms of area 
mean values of the same blocksize as the anomaly data, one can use 
the values 

Ag*tj = Agjj + tc^ (2.157) 

as "observations" in the adjustment, considering these values to refer 
to the surface of the ellipsoid. Of course the systematic reductions 
<Sgf j have to be applied in this case as well. 


CHAPTER III 


DATA USED IN THIS STUDY 


3.1 Data Requirements for Surface Gravity Analysis 

Estimation of harmonic coefficients of the earth’s gravitational 
potential based on the mathematical formulation presented in Chapter 
II, requires two types of surface observations; area-mean values of 
free-air gravity anomalies and area-mean values of orthometric heights. 
As it was mentioned in the introduction, this study aims to provide 
optimum estimates of the geopotential based on terrestrial measurements 
alone . Accordingly, altimetric observations as well as any other type 
of observations obtained from satellite techniques will not be 
considered here. 

A brief description of the fundamental data bases that were used 
in this study is given in the following sections, mainly focusing on 
some of their characteristics that are important to consider for the 
appropriate use of these data and for accurate interpretation and 
evaluation of the results. 


3.1.1 The June 1986 Gravity Data Base 

The surface free-air anomalies used in this study are taken from 
the "June 1986" data base [Despotakis, 1986], which is the most recent 
update of the OSU gravity data bank. The June 86 anomaly field 
contains 48955 l’xl* surface mean free-air anomalies based on the 
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GRS ’67 normal gravity formula [IAG, 1971]. 
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The compilation and updating of global anomaly fields, as the 
above, is based on the collection and processing of gravity information 
referring to smaller (regional) areas. Such data subsets ("data 
sources") contain gravity material in various forms, both in terms of 
data type (point or area-mean values at various subdivisions, of 
free-air or Bouger anomalies), as well as in terms of format (maps, 
tabulated data, magnetic tapes). To obtain from these sources a 
uniform global set of mean free-air anomaly estimates, some kind of 
prediction technique (the term prediction is used here in a broad 
sense) has to be applied to the data of each source. Depending on the 
nature of the source, the estimation technique varies accordingly, 
ranging from visual estimation performed on maps, to least squares 
collocation performed usually to point data or mean data of smaller 
subdivisions. The task of obtaining a reliable global anomaly field 
becomes more complicated by the fact that valuable information related 
to the original observations from which a data source has been 
compiled (such as the vertical datum to which the gravity anomalies 
refer or the exact way of computation of normal gravity), is usually 
not provided. The validity and reliability of the gravity anomaly 
estimates is deduced from comparisons with previous estimates 
(wherever possible) or "independent" data sources, such as satellite 
altimetry derived anomalies (in oceanic areas), or even from the 
experience and intuition of the person performing the compilation. 

Of equal importance with the anomaly estimate itself, is the estimate 
of its accuracy, which is deduced in most cases in a more or less 
empirical fashion. Although it is not reasonable to claim that the 
anomaly estimates in global anomaly fields are uncorrelated, no global 
field exists at present for which an error covariance function is 
defined. For regional and more homogeneous data sources, error 
covariance functions have been developed and tested [Weber and 
Wenzel, 1982]. However, the inhomogeneity of the data from which 
global anomaly fields are compiled, as well as the size of these data 
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bases introduce serious problems on the estimation of global error 
covariance functions. In this study the error covariance matrix 
associated with the June 86 anomaly field will be assumed diagonal. 

Of particular interest is a subset of 5684 values (about 12%) of the 
June 86 field, that consists of l*xl* anomaly estimates which have been 
derived by means of geophysical correlation (Figure 9). Detailed 
discussion on that estimation technique is given by Wilcox (1966, 1974, 
1978). These anomaly estimates are of inferior quality as compared to 
observed gravity (Figure 12). The poor quality of these data, in 

conjunction with their geographical distribution (they cover mostly the 
extended areas of USSR and China), creates a quite unfavorable 
environment for the estimation of long wavelength features of the 
gravity field. 

In Table 2 certain statistical quantities referring to the June 86 
field are given, while Figures 8 through 12 illustrate some aspects of 
the geographical distribution and quality of these data. 


Table 2. Statistical Characteristics of the June 1986 Anomaly Field 



Northern Hemisphere 

Southern Hemisphere 

EQ333B3 

No. of anom. 

28005 

20950 

48955 

Mean value 

-0.2 

-0.6 

-0.4 

RMS value 

28.1 

26.6 

27.5 

Min value 

-270 

-222 

-270 

Max value 

340 

192 

340 

Wtd. mean 

-0.9 

0.4 

-0.3 

Wtd. RMS 

28.5 

26.1 

27.5 

Min. std. dev. 

1.0 

1.0 

1.0 

Max. std. dev. 

62 

47 

62 

RMS std. dev. 

14.0 

14.6 

14.2 
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Figure 8. Location of the 48955 l*xl* Anomalies of the June 1986 Field. 
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Figure 9. Location of the 5684 Geophysical Anomalies of the June 1986 Field. 
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ro 10. Location of the 43271 Anomalies of the June 1986 Field that Exclude 5684 Geophysical 
Anomalies 
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In the course of this study three particular issues related to the 
gravity data have been reconsidered and examined due to their 
influence on the results; namely 

(i) The exact means of computation of the normal gravity r Q (see 

equations (2.67) and (2.69)), which is necessary to define the 

surface free-air anomaly. 

(ii) The conversion of the anomaly data from one reference system 

(e.g. GRS '67) to another (e.g. the system defined by the 

constants used in the development of the GEMT1 satellite-only 

field [Marsh et al, 1987]). 

(iii) The frequency content of the l*xl* mean anomalies of the June 
86 field. 

(i) 

In section 2.3.1 the definition of the surface free-air anomaly and 
the computational formulas required for rigorous mathematical modeling 
were given in equations (2.67) through (2.76). Before using the June 
86 anomalies it is necessary to examine whether these data comply with 
these requirements. Although for most of the gravity sources the 

relevant information is not provided, judging from the most reliable 
and accurately documented ones (e.g. Source 91 - see Despotakis, 

1986), one has to conclude that the second order vertical gradient of 
the normal gravity is neglected in the definiton of the free-air 
anomaly. Consequently, in this study it will be assumed that the June 
86 field contains area averages of anomalies defined by 


Ag = g p - [tq 0 + (fj) o «] 


(3.1) 


Probably the terms y Q and 


fill 

lahJ q 0 


have not been evaluated as 


prescribed in equations (2.69) through (2.76) (for example the term 
Nj?| is most probably assumed constant, equal to 0.3086 mgal/m) ; 

v diiJ Q q 

however these approximations will be neglected. 


The neglect of 
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the term Q H 2 though, cannot be ignored. This term can reach 

a value of about 4 mgals in regions such as the Himalayas (H ~ 7000m), 
and its effect is certainly non-neglible in view also of other more 
refined corrections (e.g. e h ). To examine its effect the following 
computations were performed: 

(a) Evaluate a global set (64800 values) of the quantity 

«gh 2 = "37 q 0 (^) 2 (3.2) 

The evaluation of 7 q 0 is performed on the mid latitude of each 
block. 

(b) Perform a harmonic analysis of the above global set, to define a 
corresponding set of harmonic coefficients. 

(c) From these harmonic coefficients the implied undulation corrections 
<5N h a can be evaluated in terms of either point or area-mean 
values. Figure 13 illustrates the corrections <5N h 2 evaluated on a 
5*x5* grid from a harmonic set corresponding to < 5 g h 2 which is 
complete up to N max = 36 (These corrections are to be added to the 
undulation obtained by omitting < 5 g h 2 ). 


LONGITUDE 



Figure 13. Undulation Correction < 5 N h 2 g asec j on a 5 « x 5 « q p1c j 
(Contour Interval is 20 cm). 
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As it can be seen, the effect reaches a maximum (absolute) of about 
-1.80 m in the Himalayas. The RMS magnitude of fig h a (up to N max = 
36) is 0.12 mgals, while the RMS <5N h a is 0.22 m. In view of its 
systematic character and its magnitude (compare with e h ), it is clear 
that such a term has to be included in the modeling. 

(ii) 

The normal equations obtained from the analysis of the surface 
gravity data are to be combined with the corresponding normal 
equations from the analysis of the satellite data. For that reason, a 
common reference system has to be adopted in the analysis of both 
data sources. In this study such a system will be defined by the 
"TOPEX constants" [Marsh et al, 1987 Section 4.1.4], i.e. 

a = 6378137. m 

1/f = 298.257 

GM = 3986004.36x10° m 3 /sec 2 
<j = 7.292115xl0~ 5 rad/sec 

To convert the l*xl* anomalies of the June 86 field, from the GRS '67 
to the TOPEX system, the following formulas are used 

Agij(TOPEx) = 4g i j( G Rs 67) + (3.4) 

i 

where 


(3.3) 


figi j = d 0 + d 2 sin 2 $j + d 4 sin 4 $ i 
and the coefficients d 0 , d 2 , d 4 are given by [Moritz, 1984] 


r ®GRS67 

r °T0PEX 

(I 2 = 3 2 

a 2 

2 GRS67 

2 topex 

d 4 = a 4 

- a 4 

GRS67 

TOPEX 


(3.5) 


(3.6) 


and 
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a 2 = ^ e 2 + k 
a 4 = | e 4 + | e 2 k 

with 



In terms of numerical values one has 


7e TO p E x = 9-7803252176 m/sec 2 

d 0 = -0.6762 mgals 

d 2 = -0.0761 mgals/rad 2 

d 4 = 0.0007 mgals/rad 4 


(3.7) 


(3.8) 


(3.9) 


Directly related to the reference system used, is the definition of the 
even degree zonals C*£ f0 of the disturbing potential. The following are 
defined for this study: 


Cfj.o = C a j , 0 + J a j(4J+l)-* i = 1,2,3 
C s£,o ~ ^a£,o £ > 4 

with [Moritz, 1984] 

I 1 -' + 5i £) 

and 


(3.10) 


(3.11) 


(3.12) 


The numerical values implied for J 2 , J 4 , J 6 using the TOPEX constants 


are 


J 2 = 0. 1082631478xl0 -2 

J 4 = -0. 237091961x10-* 
J 6 = 0. 60835 lxlO -8 


(3.13) 
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The significant digits given here for the zonal coefficients of the 
normal potential do not imply any truncation of these values. Here the 
constants given in (3.3), exact to the given digits were used to define 
the reference system, along with the above equations. The derived 
values of J a , J 4 , J 6 were transmitted to NASA/GSFC along with the 
normal equation sets. There, prior to the combination solutions the 
right side vectors of all normal equation subsets were translated in 
order to refer to a unique reference system. 

It should be mentioned here, that in the course of this study 
various other sets of constants (such as the OSU 1986 [Rapp and Cruz, 
1986a, p. 4] or the GRS ’80 [Moritz, 1984]) have been used for the 
purpose of various experiments (simulated data tests etc). Although 
these tests could have been performed equally well using the TOPEX 
constants and their results are hardly affected by the change of the 
reference values, to avoid confusion the reference values used for each 
experiment discussed in the sequel are explicitly stated. 

(hi) 

As it will be discussed in more detail in the next chapter, in the 
harmonic analysis of an incomplete set of data on the surface of the 
earth, it is of critical importance that the frequency content of the 
signal to be analysed be known, at least to an approximate degree. 
The frequency content of each l'xl' mean anomaly estimate is primarily 
a function of the number and distribution of the more detailed data 
inside the block, from which the l'xl' mean value has been estimated, 
as well as, of the estimation method (viewed as a low-pass filter) used 
to define the l'xl* mean value. To illustrate the problem, consider the 
signal being continuously known over the surface of the earth, which 
for the purpose of this discussion can be viewed as a sphere. Then it 
would be possible to design an averaging operator (low pass filter) 
with such response (transfer function), so that any desired band of 
frequencies would be eliminated from the mean values. However, in 
fact the type and distribution of detailed data inside each l'xl' block 
varies, and so does the prediction techniques used between different 
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data sources. Hence the frequency content of each area-mean value in 
a global anomaly field should not be expected to be the same. 

The consequences of that problem have been experienced in the 
current study as will be seen in Section 5.2.5. 

3.1.2 The Mean Elevation Data 

In the course of this investigation three elevation data bases have 
been examined: 

(i) OCT85.MELEV 

(ii) TUG86, and 

(iii) TUG87 

(i) OCT85.MELEV 

This digital terrain model (DTM) is a global complete set (64800 
values) of l*xl* mean elevations and has been compiled based on the 
DMA 79 model, after gross error corrections have been applied by R.H. 
Rapp. This set accompanies the June 86 gravity data [Despotakis, 
1986], since at that time it was considered the most up to date 
elevation data base. 

This set was not used in any of the computations involved in this 
study. 

(ii) TUG86 

This set was compiled by H. Sunkel at the Technical University of 
Graz and became available at OSU on July 11, 1986. It is also a 
complete set of l*xl* mean elevations. It has been compiled based on 
the OCT85.MELEV for the continental areas and the 5'x5' SYNBAPS 
bathymetric data for the oceans. This set was used in almost all of 
the preliminary investigations in this study. 

(iii) TUG87 

This set has the same origin and author as the previous one and 
became available at OSU on July 8, 1987. It has been compiled based 
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on two DTMs provided by the National Geophysical Data Center in 
Boulder, the ET0P05 and the DBDB5 (formerly called SYNBAPS). Details 
on the set up and characteristics of the TUG87 are given by Wieser 
(1987). At present TUG87 (which became available in the form of 5'x5', 
30'x30' and l*xl* mean values), is considered the most reliable source 
of elevation information. As soon as it became available, it has 
replaced TUG86 in the computations made for this study. 

In Figure 14 the locations of the 5911 l*xl* blocks where the 
absolute difference between the TUG87 and the OCT85.MELEV elevation 
exceeds 1000 m are shown. 

It should be mentioned here that the gravity anomaly and elevation 
data used, should be consistent (see equation (2.69)). However, in the 
absence of elevation information related to the gravity data sources, 
one can only use the most reliable elevation data available with 
reservations concerning their compatibility with the gravity data. 



Figure 14. Locations of the 5911 l*xl* Blocks Where IH,j TUGB7 - 
Hi joCTas.MELEV* > 1000 m. 
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CHAPTER IV 


ESTIMATION OP HARMONIC COEFFICIENTS OF THE GEOPOTENTIAL 


4.1 Combination Solutions 

The complementary spectral sensitivity of satellite versus terrestrial 
observables with respect to the gravity field, necessitates the 
combination of both data sources in order to obtain optimum estimates 
for those harmonic coefficients that correspond to frequencies 
recoverable by both techniques. The fact that satellite observations, 
such as ranges, range rates etc., can be safely assumed uncorrelated 
with terrestrial measurements, such as gravity anomalies, permits the 
combination solutions to be performed in a sequential mode. First two 
independent estimates of a potential coefficient set are obtained, one 
from the analysis of satellite data and one from surface gravity data. 
Then, based on these estimates and their associated covariance matrices 
an adjustment can be performed in order to yield a unique (with 
respect to the optimization principle adopted) set of coefficient 
estimates. 

Estimates of harmonic coefficients from orbital analyses are obtained 
from general dynamical solutions. In these solutions several subsets of 
observations (different satellites and data types) are used, and the 
potential coefficients constitute only a subset of the total set of 
parameters being estimated. This set also includes station coordinates, 
earth orientation, tidal and orbital parameters. 

The modeling and estimation techniques used is such solutions will 
not be discussed here, since these topics are beyond the scope of the 

65 
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present study. Elaborate discussions on related topics can be found in 
Kaula (1966a), Cappellari et al. (1976), Martin et al. (1980) or Marsh et 
al. (1987). For the purposes of this discussion it suffices to consider 
a simplified, schematic description of the adjustment process used in 
satellite solutions. 

Let X a denote the set of potential coefficient parameters, while Y a 
the set of all other parameters (orbit, earth orientation, tides). Then, 
linearization of the functional relations of the satellite observables with 
respect to X a and Y a leads to observation equations of the form (see 
also Appendix B for the notation) 

V S = BY + AgXg + Lg (4.1) 

where 

Lg = Lg - Lg (4.2) 

and Lg denotes the approximate values of the observables computed 
from the approximate values X° and Y° of the parameters. Matrices A s 
and B contain the partial derivatives of the observables with respect 
to X a and Y a (evaluated at X° and Y°). In principle, the spectrum of 
the gravitational field extends to infinity so that matrix A s should have 
an infinite number of columns. However, at satellite altitudes, only the 
lower degree harmonics perturb the orbits to the extent that these 
perturbations can be observed and separated from measurement noise. 
Therefore, for satellite applications, it is justifiable to truncate the 
spectrum of the gravitational field to a finite degree. 

If P 8 denotes the weight matrix associated with L| then 
minimization of the weighted norm of the residuals v|P s v s under the 
condition (4.1) leads to the normal equation system 
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b t p s b 

A^P a B 


b t p b a s 

A a P s A a 



• A 

Y 


b t p s l s 


A 



LxJ 


aIp b l s . 


(4.3) 


Elimination of the Y unknowns from the system (4.3) will result in a 
system 

N S X S = U s (4.4) 

where 

N b = [A a P a A s - (A a P a B) (B T P s B) _ 1 (B T P 8 A s ) ] (4.5) 

U s = [(A a P a B)(B T P 8 B) _1 (B T P a L a ) - A a P s L a ] (4.6) 

A modification of the above procedure was used in some GEM models 
[Lerch et al., 1979]. The modification consists of the use of a-priori 
information for the geopotential coefficients in the form of an a-priori 
weight matrix P x which is diagonal with elements equal to the reciprocal 
degree variances per coefficient as implied by a specified degree 
variance model (Kaula’s rule of 10“ 5 /n 2 ). This scheme corresponds to 
the case of "weighted parameters" [Uotila, 1986, Section 5.5] where the 
target function to be minimized is not v|P a v s but 

v I p s v s + v x P x v x = min (4.7) 

with v x denoting the residual associated with the "observed" parameters 
L x (see also Appendix B). The effect of this modification in the 
formulation given above is that the submatrix AjP a A s of the normals has 
to be replaced by (A|P s A a + P x ). The advantage of this procedure is 
that (AjP a A a + P x ) is better conditioned than A|P s A a , a fact which 
improves the estimability of certain (higher degree) coefficients which 
are highly correlated [Lerch et al, 1979]. In this particular application 
of this procedure the observed values of the parameters, L x , are assumed 

A 

equal to zero, so that v x becomes identical to X s [Marsh et al*, 1987, 
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Section 8. 1] . 

In any case, relevant to the folio-wing discussion is that a set of 
normal equations as (4.4) can be obtained from the orbital analysis. 

If a corresponding set 

N t X t = U T (4.8) 

is formed based on terrestrial gravity information, then, in principle, 
the combination solution can be obtained by 

X c = (N s + N t ) _1 (U s + U T ) (4.9) 

Of course, for (4.9) to be valid the same approximate values have to be 
used in both the satellite and terrestrial adjustments. Also, questions 
of relative weighting between N s and N T may have to be considered. 
However, these issues do not change the general principle and strategy 
of the combination solution as presented above. In the remaining part 
of this chapter various aspects related to the formation and structure of 
the normal system (4.8) are examined. 


4.2 Estimation of Potential Coefficients from Surface Gravity 

The problem of estimating a set of harmonic coefficients of the 
geopotential from surface gravity data is formulated as follows: 

Given : 

(a) A set of equi-angular area-mean values of surface free-air anomalies 
(June 1986 field), which constitutes a subset of the Nx2N (N = 7r/AX) 
set that covers the entire area of the ellipsoid. 
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(b) A diagonal S L b matrix containing estimates of the error variances of 
the above data. 

(c) A complete (Nx2N) set of mean orthometric heights. 

From these data, and based on the mathematical model (2.142) and the 
systematic reductions 5g h a and Sg r as defined in equations (3.2) and 
(3.5), one seeks to estimate a set of harmonic coefficients of the 
geopotential, complete to degree and order N max , where N max < N. For 
reasons that will be explained in the next section, it is assumed that one 
can evaluate the contribution to each data value of all harmonics of 
degree n > N max . Hence, define 

(«*hf)ij = TT H 2 ? I (n-l)M-)" l (cJ m IcJ + S n- ,Isi)IPi„ 

J a <r, r lj n=Hmax+i m=0 

(4.10) 

so that one can write (see also equation (2.142)) 


= 


1 

A<r^ 


fiM Nmax f a ) n rt 

m I (n-l)r l 

m=c 


- 2 U 

Si? 


(cj„ici 


+ Sn.isijipi, 


(4.11) 


where 

dgjj = Ag 1 j(j un0 e«) + 5g?j + (Agh 2 )ij + “ (Aghf)ij (4.12) 

Since the vertical datum inconsistencies have not been modeled in the 
solutions made for this study, the term AWfc and the superscript k in 
the above equations have been omitted. 

Equation (4.11) is identified as being of the general form 


L a = F(X a ) 
where 


(4.13) 


X“ = 


r c* r l n < 

L • • • ^nm • • • J lo( = 


^roax » 

0,1 


m < n 


(4.14) 


with 
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C* 

^nn 


if <* = 0 
if <x = 1 


(4.15) 


while the observations L b are considered to be 


1 — [ « . * dg ^ j * » ♦ ] . 


(4.16) 


The function F, a8 it can be readily seen from (4.11) is linear with 
respect to the parameters X a . Hence following the linear observation 
equations adjustment model [Uotila, 1986, Section 3.2.1], and assuming 
approximate values equal to zero, one has a set of observation equations 


v = AX a - L fc 


(4.17) 


Hence, with 


P = 4 I L b 


(4.18) 


minimization of the weighted norm of the residuals under condition (4.17) 
yields the normal equation system 


(A T PA)X* = A T PL b 


(4.19) 


and thus the least squares estimate of X a is given by 


x a = (A T PA) VPL b 


(4.20) 


while the covariance matrix of X a is given by 


I; = *o(A T PA) ‘ 


(4.21) 
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In terms of a geometric interpretation, it is well known that, if an inner 
product is defined over the column space of the design matrix A (A 
must have full column rank) by 

<x,y> p = x T Py (4.22) 

A 

then, the estimate X s above, represents the projection of L b onto the 
column space of A. 

The estimator Et B - (A T PA) _1 A T P becomes also the best linear 
unbiased estimator if the noise n of the observations has zero mean 
(E{n}=0) and S L b is defined by 

X L „ = E{nn T } (4.23) 

In addition if the probability distribution of n is Gaussian, then Ej a 
represents the maximum likelihood estimator as well [Colombo, 1981]. 


4.2.1 Structure and Characteristics of the Normal Equations 

Although the presentation given in the previous section completely 
defines the adjustment procedure used in this study, it is useful to 
re-examine some of the equations given before when these are written in 
explicit analytical form rather than in compact matrix notation. Such 
considerations are not only useful because they provide insight on the 
structure and characteristics of the normal system (4.19), but, become 
necessary also since the large size of certain matrices involved in the 
adjustment process prohibits the direct use of matrix operations. 

Analytical expressions for the elements of the matrices of the normal 
system can be derived in a compact form useful for further 
interpretation, only if the effect of data gaps (locations where no 
observation is available) is considered in an indirect manner. In the 
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following derivations an empty block is considered to contain an 
undefined data value with weight equal to zero. Such treatment of data 
gaps enables a compact and patterned formulation, efficient for 
numerical implementation, and yields identical results with the straight 
forward approach. It should be emphasized that this is to be 

interpreted only as a tool to simplify derivations and computational 
algorithms. It has no impact on the positive definiteness of the weight 
matrix of the observations or the degrees of freedom of the adjustment. 

From (4.11) it can be seen that the elements of the design matrix A 
are given by 


[A]'* 

Wnn 


1 GM 
A<r, r,j 2 




where 


v°"j - 

1 nm “ 


IPnmlC^ 

ipLisi 


if « = 0 
if a = 1. 


(4.24) 


(4.25) 


For an arbitrary element of the normal matrix N = A T PA corresponding 
to the unknowns C^ m and c£ s one has (since P is assumed diagonal) 


[N] c <* p 


N-i JN-l Ij H 

I l [A] 5 [A ]J Pij 

1=0 j = 0 ^nm L rs J 


(4.26) 


or 


[N] r ot = GM 2 (n-1) (i 


N-i 

-1) I 


IP 1 

nm 


A tf i 2 j=o 


t \ n+r 

—1 

fici 1 . 

i ci 

I 

iisi 

>is{. 


PU 
(4.27) 


For a diagonal element corresponding to C# m one has 
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[ipL 

| 2 2 y _1 

1 f a ] 


l A<f{ • 

1 J=° 

fij 4 ^ijJ 

u 


Pij (4.28) 


while the corresponding element of U = A T PL b is given by 


[U] r c( 

'-'nm 


GM(n-l) 


Y *&• ’Y r^(r-)" 

i -o A<r< j — o rij ^rjj^ 


ICjJ 

isj 


Pi jdg< j. 


(4.29) 


Consider now an idealized case where the observations are given on 
the surface of a sphere of radius r = a, the Nx2N data grid is full (i.e. 
no data gaps exist), and all observations have the same standard 
deviation <r, so that P = <r~ 2 I. Without loss of generality further assume 
<f - 1 to simplify the following expressions. In all subsequent discussion 
it is assumed that the coefficients being estimated, form a complete set 
up to degree and order N„, ax where N roax is always smaller than N = 
n/AX, the Nyquist frequency implied by the data sampling. In such case 
the previous equations become 



(4.30) 


[N] ot fl 

Kj nar'r* 


r 2 (n-l)(r-l) 


N-l 

I 

1=0 


IP 1 IP 1 
x ^njn xx rs 

Afft 2 


J=0 lisJJlisJJ 


(4.31) 


[N] r c< 


2 

y 


(n-l) a 


tzL] 

i’ *f*| 


l A<r,J 

j=« 1 

'isi 


(4.32) 


and 


. . N— i Tp 1 

[U]„« = y(n-l) l ^ 


; nm 


i =o 


2N-X 

I I 

1 j=o 


ici 

isj 


dgij 


(4.33) 


where 
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(4.34) 


However, it is true that the orthogonality of the sine and cosine 
functions is maintained in the interval [0,2rr) when discrete mean values 
of these functions are regularly sampled, provided that the Nyquist 
frequency is not exceeded. Hence, 


JN-l 


Z IC^IC^ =0 if m * s < N 


J = ° 


2N~1 i t 

l is^isJ 

J=° 


2N-1 i * 

I IC^ISJ 

j = 0 


= 0 


= 0 


if m * s < N 


for all m and all s 


(4.35) 


Equations (4.35) can be proven easily if one recognizes that [Colombo, 
1981, p. 4, equation 1.7] 


( is!)' = {-B(m)) COsmjAX + ( A(H)) siranjAX 


A(m). 


where 


A(m) 


B(m) 


sinmAA 

m 

AA 

cosmAA-1 

m 

0 


if m*0 
if m=0 
if m*0 

if m=0 


(4.36) 


(4.37) 


and makes use of the relationships [ibid, p. 10, equation 1.20-a] 
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2N-1 

1 cosmjAXcospjAX = 0 

J=° 

if m * p < N 

2N-1 


I sinmjAXsinpjAX = 0 

J=o 

if m * p < N 

2N-1 


I cosmjAXsinpjAX = 0 

J=0 

for all m and all p 


In addition the equatorial symmetry of the grid (A<r j = A<r N _ ( ) in 
conjunction with the parity properties of the integrated values of the 
associated Legendre functions, namely 

IPL = (-l) n ~ m IP„!T , ~ 1 (4.39) 

yield 


n— » TP 1 ip 1 

I — 2 r a = 0 if n-r is odd. (4.40) 

1 =o Afff 

Hence, from (4.31), (4.35) and (4.40) follows that (see also [Colombo, 1981, 
equation 2.85]) 


[N} r c< = 0 

^nrrr-'rs 


if 


<x * ft , m * s 
or 

n - r is odd 


(4.41) 


while if neither of the conditions in (4.41) apply then [N]C# m c£ s may or 
may not be zero. 


Since an element of the normal matrix N also represents the inner 
product of the corresponding columns of the design matrix A and since 
these columns consist of successions of the area averages of the surface 
spherical harmonics (see equation (4.30)), the relation (4.41) implies that: 

(a) The area averages of the surface spherical harmonics do not 
constitute an orthogonal set of base functions on the sphere, even 
if the sampling is regular and the Nyquist frequency is not 
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exceeded. The non-orthogonality (or skewness) of this set of base 
functions implies that if a function (such as the gravity anomaly) 
defined by discrete, regular, area-mean samples of it (such as the 
dg|j), is spectrally decomposed with respect to the above set of 
base functions, using quadrature formulas (see Section 4.2.2), then, 
synthesis of these spectral components will fail to recover the 
sampled values. The discrepancies between the original and the 
synthesized samples is due to the omission of the existing 
correlations between the base functions when quadrature formulas 
are used (see equation (4.47)). These correlations arise from the 
use of a finite sampling interval and occur not only in the case of 
area-mean samples (as stressed by Mainville (1986, p. 105)) but also 
in the case of regular, point value samples as it is formally shown 
by Colombo (1981, p. 53). 

(b) Since the normal matrix N is not diagonal, the estimated coefficients 
from a least squares adjustment are correlated (even if the sampling 
is regular and the Nyquist frequency is not exceeded). Hence, if a 
complete set of coefficients up to N max is estimated from a noiseless 
signal (simulated data) containing frequencies beyond N max , then the 
estimated coefficients will be distorted with respect to the "true" 
ones. This effect is called aliasing. However, it should be 

emphasized here that the cause of it is the finite block size which 
destroys the orthogonality of the surface harmonics and not the 
exceeding of the Nyquist frequency by the set of coefficients being 
estimated. In Section 5.2.4 the results from a number of 
experiments performed in order to illustrate these aspects are 
presented. 

We should note also that the conclusions drawn above are valid not 
only for the sphere, but in general for any surface of revolution 
symmetrical about the "equator" (e.g. an ellipsoid of revolution), 
provided that the data grid is symmetrical with respect to the equator 
and AX is constant. Also the use of a unit weight matrix is not 

necessary. The same conclusions could have been drawn using weights 
for the anomaly data that are longitude independent and symmetric with 
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respect to the equator, i.e. p, = p N _ wl , 

Proceeding from the idealized situation examined above to the "real 
world" case corresponding to the problem at hand, it is realized that 
the patterned structure of the normals will be destroyed once data 
gaps, unequal weights and data referring to the earth's surface are 
introduced. The dominant factor affecting the structure of the normals 
is the data gaps. Their arbitrary locations make it practically 
impossible to predict their effect on the elements of the normal system. 
If a rigorous adjustment procedure is to be used then the entire normal 
matrix (upper triangular part) has to be formed. Otherwise, if an 
approximate solution is sought, considering only a limited bandwidth of 
the normal matrix, a particular ordering of the unknowns has to be used 
so that the largest (in absolute sense) off-diagonal elements can be 
brought closest to the diagonal. Such ordering may approximately be 
determined by examining the pattern of the off-diagonal elements from 
test solutions (Wenzel, 1985]. To illustrate the effect of the ordering of 
the coefficients in conjunction with the data gaps, on the structure of 
the normal matrix, the following examples were performed. We 
considered six different ordering patterns of the unknown coefficients, 
as schematically shown on Table 3. For each particular ordering pattern 
two sets of normal equations were formed, both referring to an 
expansion up to N max s 12. In both cases the data were assumed to 
refer to the surface of the reference ellipsoid and to have the same 
accuracy (P=I). However, the first set of normals is formed from data 
that are assumed to occupy the 48955 locations of Figure 8, while the 
second from data assumed to occupy the 43271 locations of Figure 10. 
From these normal matrices the corresponding correlation matrices were 
formed. 

In Figure 15 the structure of these correlation matrices for two 
particular ordering patterns is illustrated, for each case of data 
coverage. Each "x" denotes a correlation coefficient that is larger than 
or equal (in absolute value) to 10%. The number of such values in 
cases (a) and (c) is 1569 while in cases (b) and (d) is 3724. The 
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Table 3. Ordering Patterns for the Unknown Coefficients 



ORDERING 

PATTERN 


I 

II 

III 

IV 

V 

VI 

coo 

COO 

COO 

COO 

COO 

COO 

C20 

C20 

C20 

C20 

C20 

C20 

C21 

C21 

C21 

C30 

C30 

C30 

C22 

C22 

S21 

C21 

C21 

C21 

C30 

S21 

C22 

C31 

C31 

S21 

C31 

S22 

S22 

C22 

S21 

C31 

C32 

C30 

C30 

C32 

S31 

S31 

C33 

C31 

C31 

C33 

C22 

C22 

S21 

C32 

S31 

S21 

C32 

S22 

S22 

C33 

C32 

S31 

S22 

C32 

S31 

S31 

S32 

S22 

S32 

S32 

S32 

S32 

C33 

S32 

C33 

C33 

S33 

S33 

S33 

S33 

S33 

S33 


advantage of ordering V over ordering I in case of a solution 
considering only a limited bandwidth of the normals is rather profound. 
Pattern V differs from the pattern discussed by Colombo (1981, p. 53) 
only by the fact that the parity of n-m is not taken into account in the 
ordering of C£ m . 

In our case the particular ordering used is really irrelevant since a 
rigorous solution considering the whole normal matrix will be pursued. 
For consistency purposes pattern I is used in all solutions here, which 
is the same with the ordering used in the formation of normals from 
satellite observations. 

Comparing the correlation matrices of Figure 15 by columns one 
recognizes that the number and locations of the existing data (48955) is 
more or less sufficient to assure diagonal dominance of the normals. 
However, the exclusion of the geophysically predicted data results in 
practically full normal matrices. This fact has an important consequence 
on the computation of anomaly degree variances from the coefficients 
obtained from solutions where the geophysically predicted data are 
excluded, as will be seen in Section 5.3.4. It must be emphasized that 




(a) Ordering I (b) Ordering I 

48955 Observations 43271 Observations 










the structure of the normals is affected not only by the number of data 
points excluded; the location of the resulting empty blocks is also 
important. 

Finally, it should be mentioned that the above results pertain to the 
case of data of the same accuracy, and thus reflect the influence of the 
geometry in the structure of the normals. However, in the more 
realistic case where the accuracy of the data varies, the structure of 
the normals will be also influenced by the variation of the weight from 
one measurement to another. 


4.2.2 The Least Squares Adjustment as Compared to Quadrature Formulas 

Numerical quadrature formulas (i.e. discretized couterparts of the 
orthogonality relations) have been proposed by Kaula (1966b) and 
extensively studied and used by Rapp (1969, 1977, 1981, 1984, 1986, 
1986a, b) for the estimation of potential coefficients from discrete mean 
values of gravity anomalies. A comparison of this technique to the least 
squares adjustment procedure was already made at the "early days" of 
global gravity modeling [Rapp, 1969]. It is considered appropriate here 
to briefly present some important aspects of such a comparison in a 
slightly different way. 

Consider that a complete (Nx2N) Bet of anomaly data dgjj is defined 
on the surface of a sphere of radius r = a. A set of coefficients C$ m 
which would minimize the quantity 


♦ 


(O 


M-l JH-1 


1 l 

1=o j=o 



•v Nona x n 

tJ" l (n-1) l 

1 n = o 
n*l 


i=o 


| (4.42) 


is sought. 

For an arbitrary coefficient C£ s the condition 



yields 
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The condition 


T T'y^y^ = 0 

1=0 j =0 
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(4.45) 


is imposed so that (4.44) becomes 


c <* - 
^rs " 


r(r-l) 


fN-i 2 N -1 (v 0 ** 2, | 1 N — 1 2 N — l -rJi 2 
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(4.46) 


Condition (4.45) (in integral form) is automatically fulfilled only if the 
sampling is continuous (infinitely small sampling interval). Hence, for 
surface harmonics of the same degree and order , one should expect that 
departures of their inner products from zero would decrease with 
decreasing sampling interval. 


Note that in view of (4.35) and (4.40) the condition (4.45) reduces to 


N l 2N l — pj \ i —pj* { 

I I = o, 

1=0 J =0 


if n-r even. 


(4.47) 


If further the condition that 


N-i 2 N -1 fv 0 * 1 J \ a 

l l ^ 1 = 4tt 

1=0 j=0 i 


(4.48) 


is imposed then from (4.46) one gets 
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47ry(r-l) i=0 j = o 


(4.49) 
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which is recognized to be identical to the quadrature formulas (without 
quadrature weights qi) UBed for the determination of potential 
coefficients from discrete mean values. Accordingly, quadrature 
formulas represent a least squares solution provided that: 

(i) A complete (Nx2N) set of gravity data dg { j is defined on a sphere. 

(ii) The weight function for these data is defined by p< j = cA<xt with c 
being any positive constant. 

(iii) Conditions (4.47) and (4.48) are imposed on the adjustment. 

Condition (4.47) amounts to setting all off-diagonal elements of the 
normal matrix N equal to zero, while (4.48) enforces the diagonal 
elements to take values 

N_c( = 4 tt 7 2 (n-1) 2 (4.50) 

Since these conditions are externally imposed on the system, the 
quadrature formulas will fail to recover a set of coefficients C# m from 
synthetic data dg^ that have been computed from these coefficients 
[Rapp, 1986]. 


Colombo (1981) has introduced "quadrature weights" q^ and has 
modified equation (4.49) to 


r 0 ' - 

°rs “ 


1 

47ry(r-l)qi 


N-i 2N-1_ cM , 

l I Vf^dg.j. 

<=0 j=0 J 


(4.51) 


He has also investigated numerically the definition of qj in an optimum 
manner (ibid, p. 76), so that the estimator of the harmonic coefficients 
defined in (4.51) minimizes the sum of squares of the coefficient errors. 


To illustrate these aspects the following numerical tests have been 
performed. 


We started with the Rapp (1981) potential coefficient set complete 
from N mln = 2 to N max = 36, and computed a complete set (64800 values) 
of l*xl* mean anomalies on a sphere of radius R=6378137 m, using the 
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SSYNTH program of Colombo (ibid, p. 101). Prom these synthetic data 
(zero noise) we have attempted to recover the potential coefficients 
using the following procedures: 

(1) Least squares adjustment with weight matrix P = I 

(2) Least squares adjustment with weights Pij = 

(3) By applying equation (4.46) 

(4) By applying equation (4.49), i.e. using quadrature formulas without 
optimum quadrature weights. 

In all cases the entire set of anomalies was used as input. The 
recovered coefficients from each procedure were then compared to the 
original ones. The criteria used in the comparisons were the percentage 
difference by degree and cumulatively, and the anomaly and undulation 
differences by degree and cumulatively (see also Section 5.2.3 for the 
definition of these quantities). In Table 4 the results of these 
comparisons are given. Along with the results of the four methods 
above, for certain degrees, the results of a corresponding comparison 
performed by Rapp (1986, p. 378, Table 2) where the recovery of 
coefficients was attempted using equation (4.51) with qi defined by 
Colombo (1981, p. 76) are listed under "Method 5". 


Table 4. Comparison of Input and Output Potential Coefficient Sets 



METHOD | 


(1) 

(2) 

(3) 

(i) 

(5) 1 

n 

<5N 

um 

p 

T5TM 

um 

p 

<5N 

um 

p 

6N 

um 

P 

6N 

um 

P 

2 

.00 

.00 

.00 

.00 

.00 

.00 

.10 

.00 

Til 

.20 



.10 

.00 

.01 

6 

.00 

.00 

.00 

.00 

.00 

.00 

.10 

.00 


.50 

ffij 





10 

.00 

.00 

.00 

.00 

.00 

.00 

.10 

.00 


.60 

8f!|| 

.24 

.10 

.00 

.06 

12 

.00 

.00 

.00 

.00 

.00 

.00 

.10 

.00 


.40 


.34 




18 

.00 


.00 

.00 

.00 

.00 

.10 

.00 

.16 

.50 

flfiH 

.74 




24 

.00 


.01 

.00 

.00 

.01 

.10 

.00 

.32 

.50 

8H3 





30 

.00 

.00 

.01 

.00 

.00 

.01 

.10 

.01 


.70 



.20 

.01 

.47 

36 

.00 

.00 

.01 

.00 

.00 

.01 

.10 

.01 

.47 

.70 

.04 





to 30 
to 36 

.00 

.00 

.00 

.00 

.00 

.00 

1.00 

.02 

.22 

.30 

.12 

1.07 

.80 

.02 

.20 


Note: Units of 5N are cm; units of 6g are mgals; units of P is % 
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The numerical results shown above merely confirm the conclusions that 
can be drawn from the theory, i.e. 


(i) The least squares adjustment procedure is able to recover exactly 
the original coefficients; the small percentage differences at degrees 
24, 30 and 36 are to be attributed to numerical noise. 

(ii) In the case of simulated data (zero noise) the least squares 
solution is independent of the selection of the weight matrix P. 
This is expected since in such cases all residuals become identically 
zero, so that v T Pv reaches its absolute minimum (v T Pv is a 
quadratic form), regardless of the selection of P. In case of real 
data, of course, the selection of P influences the final result (we 
do not consider here an overall scaling of P which does not affect 
the least squares solution but only its covariance matrix). 

(iii) The quadrature formulas, as expected, fail to recover the original 
coefficients. The discrepancies found in method (3) indicate the 
effect of condition (4.47), while comparing the results of (3) and 
(4), one can see the additional effect of imposing condition (4.48). 
A comparison of the results of methods (3) and (5) reveals that 
the introduction of q^ as defined in [Colombo, 1981, p. 76], almost 
compensates the effect of (4.48). However, as expected, the 
omission of the terms 


V *7 

i> i> 1 M 1 rm 
1=0 j=0 


with n-r even, implied by (4.47) cannot be offset, by the use of 
quadrature weights. Reversing the argument here, we can view 
equation (4.46) as an alternative way of defining the quadrature 
formulas. Such a definition will certainly not relieve one from the 
distortions of the spectrum caused by the enforcement of 
condition (4.47), but at least it will not introduce the additional 
distortions caused by (4.48) which is not needed for the 
development of equation (4.46). Note that, for computational 
purposes, the evaluation of the terms 
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V T 1 i£!hl 

iko ito A<r, 

appearing in (4.46) is required. Although it is certainly more 
demanding than the evaluation of q$, it only depends on the 
blocksize and the maximum degree of the expansion, so that, it can 
be performed once with its results stored and used repeatedly in 
subsequent applications. 

Finally, to obtain a quantitative feeling for the effect of 
condition (4.47), we have formed a set of normals corresponding to 
an expansion up to N max = 12 from a global set of l*xl* data 
assumed to refer to the surface of the reference ellipsoid. From 
this normal matrix the corresponding correlation matrix was then 
computed. The maximum (absolute) correlation was found to be 
48% referring to cf 0 ,o and cf 2 ,o* 

Let us now return to the problem of the combination solution. We 
assume that a complete set (Nx2N values) of gravity anomalies are 
defined on the surface of a sphere. These data now contain noise and 
are of varying accuracy. In addition a set of potential coefficient 

estimates is given, obtained from the analysis of satellite observations. 
This set is accompanied by its covariance matrix. The combination 
solution can be carried out in two ways which, following Rapp’s (1986) 
notation will be denoted as: 

Method A : The estimates of potential coefficients from terrestrial data 
are obtained from quadrature formulas. 

Method B : The estimates of potential coefficients from terrestrial data 
are obtained from a least squares adjustment. 

It is formally shown in Appendix B that in both cases the satellite 
information enters the combination solution in the same manner (compare 
equations B.13 and B.33). Hence, any differences in the final result of 
these methods has to be attributed to the different ways of obtaining 
the terrestrial estimates of the coefficients. These differences will be 
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due to: 

(i) The different minimization principle used in the least squares 
adjustment and the quadrature formulas. The minimazation principle 
in the least squares adjustment uses a weight function p^j = c/ff^j 2 
where <r f j 2 is the variance of each observation and c is any positive 
constant, while in the quadrature formulas p < j □ cA<r j . Accordingly 
quadrature formulas will yield coefficients that tend to reproduce 
the data regardless of their accuracy variations. 

(ii) The introduction of the conditions (4.47) and (4.48) in the case of 
quadrature formulas. Since these conditons are not fulfilled by 
area-mean samples of finite blocksize they will distort the spectrum 
of the field obtained from the terrestrial data analysis. 

Unfortunately, the only numerical comparison of the results obtained 
from the two combination methods, given the same input data, was 
performed by Rapp (1969) when the gravity data available were much 
fewer and less accurate as compared to the data available now. It 
would be of great interest to repeat such a comparison, based on more 
complete data and more refined modeling techniques. 



CHAPTER V 


NUMERICAL RESULTS AND COMPARISONS 


5.1 The Computational Task - Organization of the Numerical Analysis 

Based on the mathematical modeling and the estimation method 
presented in Chapters II and IV respectively, and using the data 
described in Chapter III, this study aims to accomplish the following 
computational goal: 

Form a set of normal equations corresponding to a low degree 
(< 50) harmonic expansion of the geopotential, in an optimum 

manner. 

This set is then to be combined with the corresponding normal 
equations obtained from satellite data analysis at NASA (Goddard Space 
Flight Center) for the development of an optimum global gravity model 
to be used for the TOPEX/POSEIDON mission [Marsh et al, 1987]. 

In terms of the notation used in Chapter IV, the end product of 
this work consists of the matrix N T and the vector U T appearing in 
equation (4.8). Obviously, in this analysis we have not restricted 
ourselves to the formation of the normal equations, but, for each set of 
normals formed, the adjustment of the terrestrial data was completed 
by evaluating the adjusted coefficients X T , the residuals, as well as the 
associated statistics (v T Pv, etc.). These information are essential 
for comparison purposes, in order to validate our results, identify 
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potential problems and investigate alternative ways of facing these 
problems. On the other hand, validation of our results cannot be 
considered complete unless the quality of the final result of the 
combination solution is also examined, by performing various 
comparisons (satellite orbit fits, Doppler undulation comparisons etc.). 

The computational effort involved in this study can be roughly 
divided into four parts: 

(i) Software Development 

The necessary software to perform the adjustment of the surface 
gravity data has been developed and thoroughly tested on the 
IBM-3081 at the Ohio State University. 

By far, the most demanding part of the calculation (both in terms 
of storage and time), is the formation of the matrix Nt and the vector 
Uj. For a scalar processor it becomes practically impossible to form 
the normal system for degrees of expansion greater than about 20, 
even if high speed large machines, such as the IBM-3081 are to be 
used. For the computations to become feasible the use of a vector 
processor is necessary. All the normal equation sets discussed in the 
following were formed on the CRAY X-MP/48 at the Pittsburgh 
Supercomputing Center (PSC). The software developed was optimized 
in order to take maximum advantage of vectorization on the CRAY. In 
Table 5 CPU times related to the performance of the software on the 
CRAY are given, for various degrees of expansion (N max ). For a 
solution to N max = 6 the formation of the normals on the IBM-3081 takes 
181.71 seconds. Hence, as it can be seen from Table 5, a speedup of 
about 60 was achieved by using the CRAY. It was estimated that the 
formation of normals up to N max = 50, which takes less than 2.5 CPU 
hours on the CRAY, would have taken more than 7 CPU days on the 
IBM-3081. 



Table 5. CPU Times for the CRAY X-MP/48 at PSC 
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CPU Time (seconds) | 

N ma x 

Number of 
Unknowns 

Formation of 
the Normals 

Inversion 

6 

46 

2.98 

.0043 

12 

166 

35.58 

.0677 

18 

358 

148.97 

.4265 

24 

622 

431.67 

1.8640 

36 

1366 

2118.51 

16.3983 


(ii) Network Setup 

From the previous discussion it is clear that the completion of this 
work involves three research and/or computation centers: 

(a) The Ohio State University (OSU), where the preparation of the 
terrestrial databases, the design of experiments related to the 
formation of the "terrestrial" normals and the analysis of results 
from the "terrestrial" point of view (anomaly and undulation 
comparisons etc.) are being performed. 

(b) The Pittsburgh Supercomputing Center (PSC), where the 
"terrestrial" adjustments are being carried out. 

(c) The Goddard Space Flight Center (NASA/GSFC), where the 
"satellite" normals, the combination solutions and the related 
comparisons (orbit fits etc.) are being performed. 

For that reason it was essential to establish efficient procedures for 
communication and data exchange between these centers. In Figure 16 
the topology of the network established and the activities of each 
center are schematically presented. 





Figure 16. Network Topology. 
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(iii) Preliminary Investigations 

These constitute the main part of the numerical analysis performed 
in this study. A number of numerical experiments were carried out in 
order to investigate various aspects of the problem, and to improve the 
procedures used in the normal equation formation. In Section 5.2 the 
numerical results and the conclusions drawn from these investigations 
are presented. 

(iv) Validation 


The investigation related to this study is still undergoing. Hence, 
any results to be presented, by no means should be interpreted as 
"final". Although at present a number of combination solutions have 
been performed at NASA/GSPC, which include normal equations formed 
from surface gravity (see Table 7); it would be unfair to comment on 
these results since they are still under testing. Hence, in Section 5.3 
we will attempt to draw some conclusions regarding the quality of our 
results from comparisons of the terrestrial only solutions with other 
existing geopotential models. 

All the results to be presented refer to expansions up to N max = 
36. Although a set of normals up to N max = 50 has been developed 
recently from surface gravity, it will not be included in the 
presentation since additional tests and comparisons are still to be 
performed regarding this set. 
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5.2 Preliminary Investigations 

In the following paragraphs the results of a number of experiments 
performed in this study are presented and discussed. To simplify the 
presentation, avoid unnecessary repetitions and provide an easy 
reference, we will begin with a collective description of the datasets 
involved in the computations, the setup for each experiment and the 
criteria adopted for some of the comparisons made. 

5.2.1 Description of the Datasets Used 

The basic datasets involved in the computations are listed on Table 
6 along with the notation adopted for their content. All these datasets 
contain values corresponding to a l # xl* equi-angular grid on the 
ellipsoid. A brief descritpion of these datasets is given next. 

<l)-<3> 

Synthetic area-mean gravity anomalies computed from the December 
1981 potential coefficient set [Rapp, 1981], complete from N min = 2 to 
N max = 12, 36 and 180 respectively. These values refer to the surface 
of the reference ellipsoid defined by the GRS *80 geometric constants 
[Moritz, 1984], whose reference field is defined by the dynamical 
constants of the same system. 

(4) 

Contains the surface mean free-air anomalies of the June 1986 
gravity data base (48955 values), discussed in Section 3.1.1. 

(5) -(9) 

These datasets contain the contribution to the l # xl* mean free-air 
anomaly from various spectral bands, as computed from different 
geopotential models. The names used are rather self explanatory of the 
content of each data set. The reference values used in the formation 
of all these datasets are the ones used in the development of the 
OSU86 series of potential coefficient sets [Rapp and Cruz, 1986a, p. 4]. 
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Table 6. Datasets Used in the Numerical Analysis 


No. 

Dataset Name 

Notation 

1 

DEC81.T012 

Sgi 

2 

DEC81.T036 

sg 2 

3 

DEC81.TO180 

sg 3 

4 

JTJN86 

dg 

5 

DEC81 . FROM37 . TO180 

hf i 

6 

0SU86D . FROM37 . T0180 

hf 2 

7 

0SU86F. FR0M37 . TO180 

hf 3 

8 

0SU86F. FROM181.TO360 

hf 4 

9 

0SU86F. FR0M37 . T0360 

hf 5 

10 

ELLCORR 

eiii 

11 

TERRAIN 

tc 

12 

FLAGSWITH 

fi 

13 

FLAGSWOUT 

f 2 

14 

TUG86 

e^Vi 

15 

TUG87 

e^v 2 

16 

0SU86F . FR0M37 . T0360 . TOPEX 

hf 6 

17 

ELLTERMS 

ei ^ 2 


( 10 ) 

Contains l’xl* mean values of the sum of the ellipsoidal correction 
terms, 

IEh J + IE^ J 

as defined in equations (2.111) and (2.112), computed from the OSU86F 
geopotential model complete from N mln = 2 to N max = 180. The OSU86 
reference values [ibid, p. 4] have been used for this computation. 

(ID 

Contains 466 values of l*xl* mean terrain corrections in the 
continental United States. These data were recevied from the Defence 
Mapping Agency/ Aerospace Center in March 1987. The locations of 
these data are shown in Figure 17. Their use is discussed in Section 
5.2.7. 



DFJOCN'AU FAGS IS 
OF POOR QTTMJTY 


9 



Figure 17. Location of the 466 l*xl* Mean Terrain Corrections Received 
from DMA/ AC. 

( 12 ) 

This dataset was formed for weighting purposes. To each 1‘xl* 
block a "flag" f ! is assigned as follows 

1 if | Agi * xl * osu 86F 2 ^ 36 - 4g 1 » xX * A djust 2 o | * 10 mgals 

= 


0 otherwise 


(5.1) 
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That is, the absolute difference of the l*xl* mean anomaly implied by 
the OSU86F geopotential model complete from N min = 2 to N max = 36, 
and the corresponding adjusted value of the anomaly implied by the 
terrestrial only solution ADJUST20 (see Table 7), iB compared against a 
threshold value of 10 mgals. The comparison was made only for the 
locations of the 48955 blocks of Figure 8. 

The rationale for the use of this information is discussed in 
Section 5.3.3. 

(13) 

Same as above but now the comparison is between the OSU86E and 
the ADJUST40 (see Table 7) implied anomalies. The comparison was 
performed only in the locations of the 43271 blocks of Figure 10. 

( 14) — ( 15) 

Contain the two l'xl* mean elevation data sets discussed in Section 
3.1.2. 


(16) 

Same as (9) but computed using the TOPEX constants defined in 
equation (3.3). 

(17) 

Contains l*xl* mean values of the quantity 
IEh j + IE^ J + IEp J 

(see Section 2.3.6), computed from the OSU86F set of potential 
coefficients complete from N mjn = 2 to N max = 180. The TOPEX set of 
constants was used in this computation. 

With the exception of (4) and (11), all other 
Table 6 are complete sets of 64800 values. 


datasets listed on 
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Computations of synthetic anomalies and ellipsoidal correction terms 
were made using the SSYNTH subroutine of Colombo (1981). 


5.2.2 The Experiments to be Presented 

In Table 7 some of the experiments performed for this study are 
listed, along with the relevant information concerning the setup of each 
experiment. With respect to this table the following comments need to 
be made: 

(i) Under the fourth column, the number of data used in each 
experiment takes three possible values with the following meaning: 
64800 : applies only in cases where synthetic data are used, since 
in none of our solutions we have used fill-in anomaly values. 

48955 : applies in both synthetic and real data experiments and 
implies that the data used refer to the locations where a gravity 
anomaly estimate is available in the June 1986 data base (see 
Figure 8). 

43271 : as above, but now data referring to the locations of Figure 
10 where geophysically predicted data are excluded. 

(ii) In the fifth column an entry zero implies that the adjustment was 

carried out on the ellipsoid. The case of ADJUST30 will be 

discussed separately in Section 5.2.7. 

(iii) In the sixth column an entry 1 implies that unit weight matrix is 

used. An entry 10 implies that the standard deviation of all 
anomalies was set to 10 mgals, while <t JUN 86 implies that the 
standard deviations of the June 1986 data base were used for 
the formation of the weight matrix (p<j = ff7j( juno«) )• The 
entries (20,38), (20,38,f,) and (20,38,f 2 ) will be separately 

discussed in Sections 5.2.6 and 5.3.3. 
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Finally, we should mention that in the cases where real data were used 
in the adjustment, the atmospheric correction, the second order 
vertical gradient correction Sgfc (Section 3.1.1) and the conversion of 
the gravity data from GRS’67 to the reference system used, were 
applied by the adjustment program. The only systematic reductions 
that need to be precomputed are the ellipsoidal corrections and the 
high frequency content of the anomaly data. 


5.2.3 Criteria Adopted for the Comparisons 


The analysis of the results obtained from the experiments listed 
above, involved, in most cases, intercomparisons between different sets 
of potential coefficients. For these intercomparisons a set of criteria 
established and used by Rapp (1986) has been adopted. These criteria 
are defined as follows: 


(1) The root mean square (RMS) undulation difference by degree and 
cumulatively. 

(a) By degree 

[ t 1 N 21% 

Z Z (AC| m ) (5.2) 

m=o <*=0 1 

where R is a mean earth radius 
(b) Cumulatively 


5N 



(5.3) 


(2) The root mean square anomaly difference by degree and 
cumulatively 
(a) By degree 


fig, = r(*-l)[z Z (AcfjT 

* L nri=0 <*=o J 


(5.4) 


I 
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where 7 is an average value of gravity. 

(b) Cumulatively 



(5.5) 


(3) The percentage differences of the coefficients by degree and 
cumulatively 
(a) By degree 


P * = 


2 

Z I (AC fj 

m-o &— o 


~~t i Z 2 * 

z z ccSL) 

m-o o(-o 


% 

100 


(5.6) 


(b) Cumulatively (average percentage difference) 


P = 


Nmax 

,z p 

t=SL 


l 


N -1 

"max x 


(5.7) 


Note that Pj and P depend on which of the two sets to be compared is 
used as a reference (denominator of equation (5.6)). In the following 
presentations the adopted convention is that the set used as reference 
is always listed first in the pair of coefficient sets being compared. 


Another two quantities that can be used to compare two coefficient 
sets are the average correlation by degree and cumulatively as defined 
by (Rapp, 1986). In our comparisons these quantities were only 
cursorily examined. 
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5.2.4 Recovery of Harmonic Coefficients and Aliasing Effects 

The experiments to be discussed next were designed and performed 
in order to investigate how the recovery of potential coefficients is 
affected by the frequency content of the sampled signal in conjunction 
with the completeness or sparsity of the sampling. 

The solutions performed for this purpose are ADJ020, ADJ021, 

SOLV18 and SOLV19, all referring to the recovery of a complete set of 
coefficients from N mjn = 2 to N max = 36. In all cases the anomaly data 
used were noiseless. The setup for these experiments has as follows: 

(1) ADJ020 : The signal contains the same harmonics as the ones to be 
recovered and the sampling is complete (64800 values). 

(2) ADJ021 : The signal contains the same harmonics as the ones solved 
for, but only 48955 values corresponding to the locations of the 
anomalies in the June 1986 field are sampled. 

(3) SOLV18 : The signal contains harmonics from N min = 2 to N max = 180 
and the sampling is as in (1). 

(4) SOLV19 : The signal is as in (3) and the sampling as in (2). 

The recovered coefficients from each solution were then compared to 
the original coefficients which were taken from the December 1981 
[Rapp, 1981] gravity model. These comparisons are given on Tables 8 
and 9. 

Examining the results from ADJ020 (Table 8) we verify that the least 
squares adjustment recovers exactly the coefficients, if the harmonics 
contained in the signal are the same as the harmonics estimated, and 
the sampling is complete. 

The results from the comparison of ADJ021 to the original 
coefficients should be carefully interpreted. First of all they indicate 

that if l # xl' mean anomalies are sampled at the 48955 locations 
corresponding to the locations of the June 1986 data base (Figure 8), 
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and if these samples are noiseless and contain harmonics from N min = 2 

Table 8. Comparison of the Coefficients Obtained From Solutions 
ADJ020 and ADJ021 With the Original Coefficients 
(December 1981 Field). 


DEC81 -VS- ADJ020 DEC 81 -VS- ADJ021 


DEGREE 

# OP COEFPS 

% DIPP. 

UND. DIPP. 

ANOM. DIPP. 

% DIPP. 

UND. DIPP. 

ANOM. DIFF 

2 

5 

0.00 

(METERS) 

0.000 

(MGALS) 

0.00 

0.00 

(METERS) 

0.000 

(MGALS) 

0.00 

3 

7 

0.00 

0.000 

0.00 

0.00 

0.000 

0.00 

4 

9 

0.00 

0.000 

0.00 

0.00 

0.000 

0.00 

5 

11 

0.00 

0.000 

0.00 

0.00 

0.000 

0.00 

6 

13 

0.00 

0.000 

0.00 

0.00 

0.000 

0.00 

7 

15 

0.00 

0.000 

0.00 

0.00 

0.000 

0.00 

8 

17 

0.00 

0.000 

0.00 

0.00 

0.000 

0.00 

9 

19 

0.00 

0.000 

0.00 

0.00 

0.000 

0.00 

10 

21 

0.00 

0.000 

0.00 

0.00 

0.000 

0.00 

11 

23 

0.00 

0.000 

0.00 

0.00 

0.000 

0.00 

12 

25 

0.00 

0.000 

0.00 

0.01 

0.000 

0.00 

13 

27 

0.00 

0.000 

0.00 

0.00 

0.000 

0.00 

14 

29 

0.01 

0.000 

0.00 

0.01 

0.000 

0.00 

15 

31 

0.00 

0.000 

0.00 

0.01 

0.000 

0.00 

16 

33 

0.00 

0.000 

0.00 

0.00 

0.000 

0.00 

17 

35 

0.00 

0.000 

0.00 

0.01 

0.000 

0.00 

18 

37 

0.01 

0.000 

0.00 

0.01 

0.000 

0.00 

19 

39 

0.00 

0.000 

0.00 

0.01 

0.000 

0.00 

20 

41 

0.01 

0.000 

0.00 

0.01 

0.000 

0.00 

21 

43 

0.01 

0.000 

0.00 

0.01 

0.000 

0.00 

22 

45 

0.00 

0.000 

0.00 

0.01 

0.000 

0.00 

23 

47 

0.01 

0.000 

0.00 

0.01 

0.000 

0.00 

24 

49 

0.01 

0.000 

0.00 

0.01 

0.000 

0.00 

25 

51 

0.01 

0.000 

0.00 

0.01 

0.000 

0.00 

26 

53 

0.01 

0.000 

0.00 

0.01 

0.000 

0.00 

27 

55 

0.01 

0.000 

0.00 

0.01 

0.000 

0.00 

28 

57 

0.01 

0.000 

0.00 

0.01 

0.000 

0.00 

29 

59 

0.01 

0.000 

0.00 

0.01 

0.000 

0.00 

30 

61 

0.01 

0.000 

0.00 

0.01 

0.000 

0.00 

31 

63 

0.01 

0.000 

0.00 

0.01 

0.000 

0.00 

32 

65 

0.01 

0.000 

0.00 

0.01 

0.000 

0.00 

33 

67 

0.01 

0.000 

0.00 

0.01 

0.000 

0.00 

34 

69 

0.01 

0.000 

0.00 

0.01 

0.000 

0.00 

35 

71 

0.01 

0.000 

0.00 

0.01 

0.000 

0.00 

36 

73 

0.01 

0.000 

0.00 

0.01 

0.000 

0.00 

AVERAGE % DIPP. 

RMS UNDULATION DIPP. 
RMS ANOMALY DIPP. 

0.00 

0.00 

0.00 

0.01 

0.00 

0.00 


to N max = 36, then these samples are sufficient for an almost perfect 
recovery of these coefficients. It should be emphasized here that this 
result pertains only for the particular blocksize used here (l’xl*), the 
particular number and location of the data gaps and the particular 
degree of expansion. If any of these "components" of the system 

changes, one should expect changes in the agreement of the recovered 
coefficients to the original ones as well, and this is due to the fact 
that the effect of data gaps is interrelated with the above 
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"components". For example, if we keep the blocksize and the location 
and number of data gaps constant, and we start increasing the 


Table 9. Comparison of the Coefficients Obtained From Solutions 
SOLV18 and S0LV19 With the Original Coefficients 
(December 1981 Field). 

DBC81 -VS- SOLV18 DEC81 -VS- SOLV19 


DEGREE # OP COEFFS 


2 

5 

3 

7 

4 

9 

5 

11 

6 

13 

7 

15 

8 

17 

9 

19 

10 

21 

11 

23 

12 

25 

13 

27 

14 

29 

15 

31 

16 

33 

17 

35 

18 

37 

19 

39 

20 

41 

21 

43 

22 

45 

23 

47 

24 

49 

25 

51 

26 

53 

27 

55 

28 

57 

29 

59 

30 

61 

31 

63 

32 

65 

33 

67 

34 

69 

35 

71 

36 

73 


AVERAGE * DIFF. 

RMS UNDULATION DIFF. 
RMS ANOMALY DIFF. 


DIPP. 

UND. DIPP 


(METERS) 

0.14 

0.026 

0.12 

0.023 

0.12 

0.012 

0.16 

0.012 

0.15 

0.009 

0.19 

0.009 

0.23 

0.007 

0.28 

0.008 

0.29 

0.006 

0.41 

0.007 

0.58 

0.006 

0.44 

0.006 

0.68 

0.006 

0.75 

0.006 

0.62 

0.006 

0.78 

0.006 

0.86 

0.006 

0.93 

0.006 

1.30 

0.006 

1.16 

0.006 

1.09 

0.007 

1.33 

0.006 

1.73 

0.007 

1.47 

0.006 

2.34 

0.008 

2.22 

0.007 

2.53 

0.009 

*.56 

0.008 

3.21 

0.011 

3.8: 

0.011 

5.79 

0.015 

5.33 

0.015 

6.91 

0.023 

29.88 

0.087 

45.12 

0.111 

3.59 

0.15 


ON. DIPP. 

% DIFF. 

(MGALS) 


0.00 

14.08 

0.01 

11.17 

0.01 

18.76 

0.01 

17.49 

0.01 

20.07 

0.01 

18.74 

0.01 

31.44 

0.01 

30.19 

0.01 

36.95 

0.01 

42.23 

0.01 

59.85 

0.01 

43.62 

0.01 

63.69 

0.01 

64.68 

0.01 

49.47 

0.01 

69.33 

0.02 

76.28 

0.02 

82.63 

0.02 

108.05 

0.02 

89.50 

0.02 

77.82 

0.02 

86.32 

0.03 

92.87 

0.02 

84.43 

0.03 

106.05 

0.03 

109.34 

0.04 

87.32 

0.04 

104.11 

0.05 

81.21 

0.05 

114.96 

0.07 

98.85 

0.08 

94.22 

0.11 

65.42 

0.45 

76.60 

0.59 

67.45 


65.58 


0.77 


UND. DIFF. 

AN0M. DIFF 

(METERS) 

(MGALS) 

2.522 

0.39 

2.112 

0.65 

1.809 

0.83 

1.294 

0.80 

1.148 

0.88 

0.911 

0.84 

0.970 

1.04 

0.817 

1.01 

0.834 

1.15 

0.719 

1.11 

0.632 

1.07 

0.633 

1.17 

0.550 

1.10 

0.522 

1.12 

0.469 

1.08 

0.524 

1.29 

0.531 

1.39 

0.516 

1.43 

0.511 

1.49 

0.455 

1.40 

0.465 

1.50 

0.395 

1.34 

0.381 

1.35 

0.369 

1.36 

0.360 

1.38 

0.348 

1.39 

0.317 

1.32 

0.333 

1.44 

0.285 

1.27 

0.319 

1.47 

0.260 

1.24 

0.272 

1.34 

0.214 

1.09 

0.222 

1.16 

0.165 

0.89 

5.02 



7.06 


maximum degree of the harmonics we estimate, (increasing consistently 
the frequency content of the sampled data), we would experience a 
degrading in the agreement between the recovered and the original 
coefficients as the maximum degree increases. In fact, at some point 
the normal system will become singular, even though the number of 
observations would still exceed the number of unknowns. The 
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particular degree at which the singularity will occur depends on the 
number and locations of the empty blocks. 

Such behavior will not be seen if the data grid is full. In that 
case we will always be able to recover the coefficients exactly (inside 
the limits of the numerical noise), as long as the maximum degree (and 
order) of the complete set of estimated coefficients, does not reach or 
exceed the Nyquist frequency implied by the data sampling. If the 
maximum degree becomes equal to the Nyquist frequency the normal 
system becomes singular, although again, the number of observations is 
larger than the number of unknowns. This is formally proved by 
Colombo (1981, p. 11). 

The results from the comparisons of SOLV18 and S0LV19 with the 
original coefficients, given in Table 9, serve as an example of the 
aliasing effect discussed in Section 4.2.1. 

In S0LV18, although the Nyquist frequency is not exceeded by the 
estimated coefficients, and although the sampling is regular and 
complete, the recovered coefficients are distorted with respect to the 
original ones. This is due to the fact that the signal has power 
beyond N max = 36 and the correlations between the coefficients from 2 
to 36 with the coefficients from 37 to 180 are neglected in S0LV18 (as 
well as in SOLV19). Schematically we can represent the normal system 
as follows 


f N n 1 

Nia 



2 

' 

1 



X, 


Ui 

2^3 6 1 




3 6 


1 

1 

N 2 2 



3 7 
4r 


1 

374180 . 


. x 2 . 

180 

. u 2 . 


(5.8) 


The exact recovery of the harmonics from 2 to 36 would have been 
given by [Uotila, 1986, p. 155] 
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X x = (N„ - N 12 N 22 nJ 2 ) l (Ux - N l2 N 22 U 2 ) (5.9) 

while S0LV18 (as well as SOLV19) represent the solution 

X? = NllUx (5.10) 

Since N i2 is a non-zero matrix, Xtf will be different from the original 
coefficients represented by X t . However, in SOLV18 the only non-zero 
elements in N 12 are the ones corresponding to C£ m and C°? m with n-r 
even (see equation (4.41)). It is thus expected that the distortion of 
the spectrum from 2 to 36 will be increasing with the degree and in 
general will be small. Severe distortion occurs close to the "cut-off 
boundary", as expected. In contrast the introduction of data gaps in 
S0LV19 masks the diagonal dominance of the normals. The 
cross-covariance matrix N 12 diverges far from the zero matrix (see also 
Figures 15a, 15c), so that its omission in (5.10) causes severe 

distortions of the spectrum. The particular degrees affected most, are 
now a function of the number and location of data gaps and the 
distortion of the spectrum does not necessarily increase with 
increasing degree any more, as it is verified from Table 9. 


5.2.5 Removal of High Frequency Content of the Gravity Anomalies 

From the results presented in the previous section it becomes 
clear, that in order to obtain estimates of the harmonic coefficients 
from an incomplete set of anomalies, without suffering large distortions 
of the spectrum due to aliasing effects, the contribution to the anomaly 
data from all frequencies beyond the ones estimated has to be 
removed. 

In principle, to do so we need to devise an ideal low-pass filter 
that will effectively eliminate all the high frequency content of the 
data without distorting the spectrum of the lower frequencies. In 
practice, to construct such a filter is impossible since the exact 
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frequency content of each l*xl* mean value of the gravity anomaly is 
unknown. There are two alternatives that can be used to approximate 
such a process: 

(1) Based on the l*xl* mean free-air anomalies, predict mean anomaly 
values over larger blocks (e.g. for expansions up to N max = 36, use 
5*x5* blocks), in an attempt to reduce the high frequency content 
of the data. 

(2) Use an existing geopotential model (which extends beyond the 

maximum degree and order for which the current solution is 

attempted), in order to estimate the contribution of the higher 

frequencies to each l*xl* mean anomaly, which then can be 
subtracted from each data value. 

Method (1) has the disadvantage that a covariance function between 
the l*xl* mean values and the larger blocksize (5*x5*) mean values has 
to be estimated, for global application, in such a way that it eliminates 

the power from all frequencies above the ones estimated, and causes 

the minimum distortion of the lower part of the spectrum. The 
difficulties involved in such a task are rather obvious. However, it 
should be mentioned, that if such a goal is achieved, the computational 
effort to form the normal system would be significantly reduced since 
much fewer "observations" would be used in the adjustment (for 5*x5* 
blocks, 2592 values provide global coverage). 

In contrast, method (2) leaves the computational effort for the 
normal equation formation unchanged, but provides much better control 
over the frequencies which are eliminated and is in general "cleaner" 
in terms of additonal assumptions and approximations with respect to 
method (1). For these reasons we have decided to proceed using 

method (2) in this study. However, before using this procedure we 
had to investigate and clarify some important aspects. These are: 

(a) The sensitivity of the results of the least squares adjustment, with 
respect to the geopotential model used for the removal of the high 
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frequency content of the data. 

(b) The maximum degree (and order) up to which the high frequency 
contribution should be computed. 

To numerically investigate the first problem we have performed three 
solutions in each one of which a different high degree field has been 
used for the removal of the high frequency contribution to the 
anomalies. In all cases this contribution was computed from harmonic 
sets complete from N min = 37 to N max = 180. These solutions have as 
follows: 

SQLV04 : Removal using the December 1981 field [Rapp, 1981]. 

SQLV03 : Removal using the OSU86D field [Rapp and Cruz, 1986a]. 
ADJUST02 : Removal using the OSU86F field [Rapp and Cruz, 1986b]. 

The intercomparisons of the results obtained form these solutions are 
given on Table 10. Before we comment on these results let us first list 
two properties that the high degree expansion is desirable to have for 
the purpose that we want to use it: 

(i) The surface gravity data used for the high degree expansion 
should be as consistent as possible with the surface gravity data 
used in the current solutions. Ideally we would like the high 
degree expansion to have been developed based on the same 
surface data with the ones used here. Unfortunately most of the 
existing high degree fields (e.g. OSU86 series of high degree 
expansions), have used, in the oceanic regions, primarily altimeter 
derived anomalies rather than ship-track estimates which are used 
in this study. 

(ii) The high degree expansion should contain a minimum of satellite 
information so that we do not introduce significant correlations 
between our reduced observables and the observables used in the 
orbital perturbation analysis. For this reason it is preferable not 
to use high degree fields where a large number of satellite 
derived fill-in values has been used in their development (e.g. 
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OSU86C or OSU86E). 

Returning now to the results of Table 10 we observe that both 
solutions SOLV03 and ADJUST02 have a substantial and about equal 
difference from the results of SOLV04. This is expected because 
OSU86D and OSU86F are almost identical in the range of degrees 37 to 
180 (see also [Rapp and Cruz, 1986b, p. 12, Figure 4]), since they have 
been developed from consistent surface gravity data, while the 
December 1981 field has been developed from a substantially different 
set of surface anomalies both in terms of number of data and in terms 
of quality. 


Table 10. Intercomparison of Terrestrial Only Solutions Where 
Different High Degree Models Were Used for the Removal 
of the High Frequency Content of the Gravity Anomalies 


DEGREE # or COEFPS 


2 

5 

3 

7 

4 

9 

5 

11 

6 

13 

7 

15 

a 

17 

9 

19 

10 

21 

11 

23 

12 

25 

13 

27 

14 

29 

15 

31 

16 

33 

17 

35 

18 

37 

19 

39 

20 

41 

21 

43 

22 

45 

23 

47 

24 

49 

25 

51 

26 

53 

27 

55 

28 

57 

29 

59 

30 

61 

31 

63 

32 

65 

33 

67 

34 

69 

35 

71 

36 

73 


AVERAGE % DIPT. 

RMS UNDULATION DIPT. 
RMS ANOMALY DIPT. 



SOLV04 -VS- SOLV03 

% DIFF. 

UND. DIFF. 

ANOM. DIFF 


(METERS) 

(MGALS) 

2.10 

0.385 

0.06 

4.40 

0.867 

0.27 

4.97 

0.523 

0.24 

7.92 

0.532 

0.33 

#.77 

0.410 

0.32 

8.66 

0.414 

0.38 

12. 11 

0.375 

0.40 

10.77 

0.309 

0.38 

15.50 

0.372 

0.51 

15.45 

0.344 

0.53 

27.34 

0.359 

0.61 

18.17 

0.296 

0.55 

29.20 

0.285 

0.57 

20.17 

0.217 

0.47 

15.07 

0.174 

0.40 

19.62 

0.203 

0.50 

24.73 

0.208 

0.54 

24.94 

0.195 

0.54 

30.46 

0.185 

0.54 

24.34 

0.180 

0.55 

23.16 

0.166 

0.54 

23.17 

0.156 

0.53 

26.48 

0.158 

0.56 

26.29 

0.149 

0.55 

27.18 

0.135 

0.52 

29.88 

0.124 

0.50 

26.77 

0.123 

0,51 

31.07 

0.130 

0.56 

26.68 

0.121 

0.54 

29.30 

0.127 

0.59 

29.30 

0.107 

0.51 

28.44 

0.106 

0.52 

17.83 

0.076 

0.38 

22.37 

0.094 

0.49 

18.18 

0.061 

0.33 

20.28 




30LV04 -VS- ADJUST02 


% DXPP. 

UND. DIPP. 

ANOM. DIFF 


(METERS) 

(MQALS) 

3.85 

0.707 

0.11 

4.02 

0.793 

0.24 

5.29 

0.558 

0.26 

7.74 

0.520 

0.32 

7.81 

0.412 

0.32 

8.24 

0.394 

0.36 

12.70 

0.393 

0.42 

11.48 

0.329 

0.41 

16.43 

0.394 

0.55 

16.14 

0.359 

0.55 

28.48 

0.374 

0.63 

19.20 

0.313 

0.58 

29.52 

0.288 

0.58 

21.08 

0.226 

0.49 

17.06 

0.197 

0.46 

21.18 

0.219 

0.54 

25.69 

0.216 

0.57 

2" .25 

0.213 

0.59 

34.06 

0.207 

0.60 

26 .08 

0.193 

0.59 

25.37 

0.182 

0.59 

24.38 

0.164 

0.55 

27.33 

0.163 

0.58 

27.02 

0.154 

0.57 

28.01 

0.139 

0.53 

30.15 

0.125 

0.50 

27.63 

0.127 

0.53 

30.45 

0.127 

0.55 

26.01 

0.118 

0.53 

29.88 

0.130 

0.60 

29.51 

0.108 

0.51 

28.43 

0.106 

0.52 

18.59 

0.079 

0.40 

22.41 

0.094 

0.49 

16.59 

0.062 

0.33 

21.06 




1.76 


2.84 


1.86 


2. 95 
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Comparison between ADJUST02 and SOLV03 yields an overall 
percentage difference of 3.52%, while the RMS undulation and anomaly 
differences are 0.58 m and 0.46 mgals respectively. However, part of 
this discrepancy is due to the fact that in ADJUST02 ellipsoidal 
correction terms have been included, in contrast to SOLV03. 

Prom the point of view of the desirable properties listed before, 
either one of the two fields OSU86D and OSU86F can be equally well 
used for the removal of the high frequency content of the data. In 
the following solutions we have used OSU86P for this removal. 

The decision regarding the maximum degree up to which the 
removal should take place was greatly facilitated by an examination of 
the residuals obtained from the adjustment of the surface anomalies. 
In Figure 18 the geographical distribution of the 11404 residuals, 
obtained from ADJUST02, which exceed (in absolute value) 10 mgals is 
shown. 

Examination of this figure indicates that the residuals are highly 
correlated with the high frequency content of the anomaly field 
(observe their signature in the area of the Andes). This observation 
leads to the following two conclusions: 

(i) The frequency content of the l'xl* mean values of the June 1986 
field extends beyond degree 180. 

(ii) The neglected correlations between the coefficients from degrees 2 
to 36 and the coefficients beyond 181 (see also equation (5.9)) are 
so small that the contribution of these coefficients to the anomaly 
data does not contaminate (at least extensively) our solution. 
Rather this contribution manifests itself in the form of residuals 
of the adjustment. 

Based on the above we have decided in the following solutions to 
remove from the data the high frequency part corresponding to the 
harmonics from 37 to 360. This part has been evaluated from the 
OSU86F set of potential coefficients [Rapp and Cruz, 1986b]. 
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A solution of this type is ADJUST20. In Figure 19 the locations of the 
9416 residuals from ADJUST20 whose values exceed (in absolute value) 
10 mgals are shown. It is interesting to note that in the land regions 
where the data used in ADJUST20 are consistent with the data in 
OSU86F, the least squares adjustment yields a remarkable fit. 

On the other hand the comparison between ADJUST02 and ADJUST20 
yields an overall percentage difference of 8.7% with the RMS undulation 
and anomaly differences being 0.74 m and 1.15 mgals respectively. 
When the two solutions ADJUST02 and ADJUST20 are compared to a 
common standard (the OSU86F field), they yield practically the same 
differences as it can be seen from the results given on Table 11. 
These numbers substantiate the argument made in (ii) above. 

An additional comment should be made regarding the computation of 
the contribution of the higher frequencies to the anomaly data. As it 
can be seen from equation (4.10) this contribution has to be evaluated 
at the surface level (r = r \ j ) in order to be consistent with the 
anomaly data which refer to the surface of the earth. However, the 
elevation introduces a long wavelength effect which can be neglected 
here since we are interested in the high frequency part of the 
spectrum. Accordingly, in this study, the terms (6g hf ) { j were 
evaluated on the surface of the ellipsoid (r = r E1 ). 

5.2.6 Weighting of the Surface Gravity Data 

There are two basic problems associated with the estimation of the 
error properties of the surface gravity data and their influence on the 
results of combination solutions: 

(a) The assumption that the errors of the surface gravity data are 
uncorrelated 
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Table 11. Comparison of Solutions ADJUST02 and ADJUST20 to the 
OSU86P Geopotential Model 


OSU86P -VS- ADJUST02 OSU86P -VS- ADJUST20 


DEGREE 

# OP COEFFS 

% DIPP. 

UND, DIPF. 

ANOM. DIFP. 

% DIPP. 

UND. DIPP. 

ANON. DIFP 

2 

5 

21.26 

(METERS) 

3.808 

(MGAL3) 

0.59 

20.18 

(METERS) 

3.615 

(MGAL3) 

0.56 

3 

7 

12.86 

2.434 

0.75 

13.60 

2.573 

0.79 

4 

9 

19.98 

1.926 

0.69 

20.41 

1.970 

0.91 

5 

11 

18.23 

1.349 

0.83 

17.92 

1.327 

0.82 

6 

13 

19.25 

1.100 

0.85 

19.05 

1.088 

0.84 

7 

IS 

21.05 

1.008 

0.93 

20.41 

0.978 

0.90 

8 

17 

35.00 

1.077 

1.16 

34.22 

1.054 

1.13 

9 

19 

31.21 

0.841 

1.03 

28.81 

0.776 

0.95 

10 

21 

34.04 

0.762 

1.06 

33.49 

0.750 

1.04 

11 

23 

51.33 

0.871 

1.34 

50.33 

0.854 

1.31 

12 

25 

45.46 

0.515 

0.87 

45.06 

0.511 

0.86 

13 

27 

37.97 

0.548 

1.01 

36.85 

0.532 

0.98 

14 

29 

63.66 

0.569 

1.14 

65.76 

0.588 

1.16 

15 

31 

63.65 

0.480 

1.03 

64.40 

0.486 

1.05 

18 

33 

71.08 

0.572 

1.32 

70.36 

0.566 

1.31 

17 

35 

78.46 

0.537 

1.32 

81.36 

0.556 

1.37 

16 

37 

66.41 

0.465 

1.22 

68.78 

0.482 

1.26 

19 

39 

71.95 

0.449 

1.24 

76.81 

0.479 

1.33 

20 

41 

62.13 

0.392 

1.14 

81.40 

0.368 

1.13 

21 

43 

91.39 

0.422 

1.30 

92.92 

0.430 

1.32 

22 

45 

66.60 

0.372 

1.20 

67.75 

0.368 

1.19 

23 

47 

70.15 

0.315 

1.07 

73.55 

0.330 

1.12 

24 

49 

77.24 

0.301 

1.06 

74.61 

0.290 

1.03 

25 

51 

61.00 

0.296 

1.10 

61.67 

0.302 

1.11 

28 

53 

79.36 

0.266 

1.02 

78.60 

0.264 

1.02 

27 

55 

56.09 

0.186 

0.74 

56.60 

0.188 

0.75 

28 

57 

51.13 

0.198 

0.82 

49.66 

0.192 

0.80 

29 

59 

66.49 

0.227 

0.98 

70.15 

0.233 

1.00 

30 

81 

50.76 

0.196 

0.88 

53.19 

0.206 

0.92 

31 

83 

53.34 

0.177 

0.81 

5 r .00 

0.182 

0.64 

32 

85 

42.41 

0.140 

0.67 

43.79 

0.144 

0.69 

33 

87 

46.34 

0.160 

0.79 

48.00 

0.166 

0.82 

34 

69 

38.21 

0.143 

0.73 

40.37 

0.151 

0.77 

35 

71 

38.94 

0.135 

0.71 

39.13 

0.136 

0.71 

38 

73 

37.13 

0.111 

0.60 

38.75 

0.116 

0.62 

AVERAGE % DIPP. 

RMS UWJULATION DIPP. 
RMS ANOMALY DIPP. 

50.73 

5.89 

5.91 

51.23 

5.82 

5.95 


(b) The need for an appropriate relative weighting scheme, so that an 
optimum combination of the satellite and terrestrial normal 
equations can be achieved. 

It was mentioned before (Section 3.1.1) that the assumption of 
uncorrelated noise for the surface anomalies is an unrealistic one. We 
recognize that especially for data belonging to the same source common 
systematic errors in the observations or the compilation procedures 
used for their estimation may be present, causing error correlations 
between such data. Weber and Wenzel (1982) have investigated 
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estimation of error covariance functions for a limited area in Europe 
based on the existing differences between dual or multiple data sources 
referring to this area. 

However, implementation of the same technique on a global basis is 
not possible at present, due to the lack of redundant data and the 
inhomogeneity of the quantity and quality of available data from one 
area to another. 

In an attempt to compensate for the omission of the existing 
correlations between the surface anomaly estimates Rapp and Cruz 
(1986a) have rescaled (increased) the variances of the anomaly 
estimates provided in the June 1986 data base, and constrained their 
values in a specified range, in order to obtain moderate weight ratios. 
The particular rescaling and the range of the variances selected have 
been determined in an empirical fashion heavily dependant on numerical 
experimentation [Rapp and Cruz, 1986a, Sections 6.1, 6.3]. 

Numerical experimentation is also employed for the determination of 
appropriate relative weights for the terrestrial and the satellite 
normals in combination solutions. The goal here is mainly to obtain 
realistic error degree variances for the potential coefficient estimates 
of the combined solution. To investigate the implications of different 
relative weighting schemes and facilitate the determination of optimum 
relative weights, a number of comparisons (potential coefficient error 
degree variances, RMS from satellite orbit fits, Doppler derived 
undulations) has been performed and examined in the development of 
the OSU86C/D geopotential models [Rapp and Cruz, 1986a, Section 6.3]. 
Corresponding comparisons provide also the means of determining 
optimal relative weights in the combination solutions performed by 
NASA/GSFC (see also [Marsh et al, 1987, p. 193, Table 8.1]). 

In this study we have performed three test solutions (only from 
surface gravity data), in order to investigate the effect of different 
weighting schemes. 
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SOLVOl : all anomaly standard deviations are set to 10 mgals 
SQLV02 : the anomaly standard deviations of the June 1986 field 
(*i j(junb«)) are used. 

SQLV03 : the standard deviation <r j j of the data value dg ( j is defined as 
follows 


j 


max (20, 2.5<r, j ( JUNa6 ) ) 
min(2.5<r, j ( JUN86 ) , 38) 


(5.11) 


Such a weighting scheme was used in the development of the OSU86C/D 
gravity models [Rapp and Cruz, 1986a]. 

The weights in all three solutions are assigned by 

Pi j = (5.12) 

where <r£ denotes the a priori variance of unit weight and has always 
been set equal to 1 (unitless). 

On Table 12 the results from the comparisons of these solutions 
with the OSU86F potential coefficient set are given. These results 
indicate no significant change on the solutions due to the various 
weights used. 

The same behavior was seen by examining the classification in 
terms of magnitude and the geographical distribution of the residuals 
from these solutions. It is interesting to note the values of the a 
posteriori variance of unit weight ($3) resulting from these solutions. 
We have 

SOLVOl: al = 1.22 

SOLV02 : = 2.51 

SOLV03 : = 0.15 


Table 12. Comparison of Solutions SOLVOl, SOLV02 and SOLV03 to the OSU86F Geopotential Model. 
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Prom the point of view of internal consistency the 10 mgal standard 
deviation for the gravity data appears to yield the best results, while 
the June 1986 and (20,38) weighting schemes indicate optimistic and 
pessimistic estimates for the accuracies respectively. 

The similarity of the results from these solutions regardless of the 
weights used should be expected. The effects of different weighting 
schemes for the anomalies will become apparent only when the dual 

satellite information is introduced in the combined solution. 

However, in our case we did not have the ability to perform 
combined solutions and examine the effect of different relative weighting 
schemes. Hence, we decided to relay on the results of Rapp and Cruz 
(1986a) and adopt (at least in the beginning) the (20,38) weighting 

scheme used in the OSU86C/D solutions. The fact that such a selection 

was not an inappropriate one was verified later, when the first 

preliminary results from combined solutions were made available to us 
[E. Pavlis, private communication]. According to these results the 
"calibration scale factors" (see [Marsh et al, 1987, p. 264 for definitions) 
by degree and by order obtained for PGS3226 (see Table 7) were 0.97 
and 0.93 respectively with their ideal values being 1.00. 

Prom the discussion given above it becomes apparent that the 

treatment of the error properties of the data is currently based on 

empirical methods heavily based on numerical experimentation. Although 
numerical experimentation (to some extent) will probably remain 
necessary as a tool guiding the analysis and the making of decisions, it 
would be highly desirable to further investigate, also from the 

theoretical point of view, the estimation of error properties for global 

surface anomaly data bases. 
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5.2.7 The Use of Terrain Corrections 

In order to numerically investigate the influence of the use of 
terrain corrections on our results we have performed the solution 
ADJUST30 with the following setup: 

ADJUST30 : the ''observations" which are input to the adjustment are 
defined by (see Table 7) 

d<j = dg{ j — (hf 5 ) j j — + tcjj (5.12) 

for the locations corresponding to the 466 blocks of Figure 17, and as 

| d^ = dg^ - (hf 5 ),j - (eiO, j (5.13) 

for the rest of the 48955 blocks of Figure 8. In the cases 
corresponding to (5.12), dg<j is assumed to refer to the surface of the 
reference ellipsoid (elevation is set to zero), in agreement with the 
formulation presented in Section 2.4. 

With the exception of the use of terrain corrections this solution is 
identical in terms of setup with ADJUST20 (weighting scheme, high 
frequency content removal etc), as it can be seen from Table 7. Hence, 
the effect of the terrain corrections is examined by an intercomparison 
of the results from these two solutions. 

The average percentage difference of ADJUST30 as compared to 
ADJUST20 was found to be 1.74%, while the RMS undulation and anomaly 
differences between the two solutions were 0.51 m and 0.27 mgals 
respectively. However, these figures are representative of the 
differences between the two fields in a global sense. Of greater 

interest here are the differences between the two solutions in the area 
where the terrain corrections were applied. 
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Figure 20 illustrates the undulation differences between the two 
solutions (N AD just2o - N A0JUST3O ) in the region of interest (♦ = 30* to ♦ 
= 50* and X = -130* to X = -70*). 


LONGITUDE 

130 -100 -70 



Figure 20. Undulation Differences - ADJUST20 Minus ADJUST30 (Contour 
Interval is 1 m, Based on a 2*x2* Grid) Superimposed are the 
Locations of Selected Laser Stations. 


Comparing Figure 20 with Figure 17 it becomes clear that the terrain 
corrections have left a long wavelength signature on the field in the 
area where they were applied. The fact that such a signature is an 
erroneous one was verified through comparisons of the undulations 
computed from the two solutions with the undulations at 13 laser 
stations in the above area. The latter, denoted by N c , were given by 
Despotakis (1987, p. 46, Table 10), and have been computed by 
differencing the ellipsoidal and orthometric heights at these stations. 
They refer to an ellipsoid with parameters 

a = 6378136.0 m ) 


1/f = 298.257222101 


(5.14) 
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The locations of these stations are also, shown in Figure 20. For the 
purpose of this comparison we have formed two geopotential models both 
up to N max = 360 as follows: 

model F20: ADJUST20 2 -» 36 augmented by 0SU86F 37 -» 36O 
model F30: ADJUST30 2 -> 36 augmented by OSU86F 37 -* 360 

The undulations implied by F20 and F30 (denoted N 20 and N 30 
respectively) at the locations of the 13 laser stations were computed and 
compared with N c . The results of this comparison are given on Table 
13. Comparing the differences given in columns five and six we verify 
that ADJUST20 is in good agreement with the "ground truth" while 
ADJUST30 suffers from the long wavelength errors shown in Figure 20. 
The last column serves as a double check of the values shown in Figure 
20 . 


At this point we have no answer to what causes the problem 
identified before. Possible causes for this problem may be: 

(a) The omission of the terms (2.148) and (2.152) may not be permissible 
for the degree of expansion we use here (36). In [Moritz, 1966, p. 
105] it is postulated that such omission is justifiable up to degree 
about 5. In our application such a guideline is certainly being 
violated. 

(b) The use of the terrain corrections only in a limited area. 

Further investigation is required in order to elucidate the appropriate 
use of terrain corrections in global geopotential modeling. 

5.3 Validation 

During the course of this study five major sets of normal equations 
have been developed from surface gravity data, to be used in 
combination solutions at NASA/GSFC. These normal sets (and the 
corresponding terrestrial only solutions) are denoted 
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Table 13. Undulation Comparisons at Laser Stations in the Western 
United States 


ID 

n g * 

N 2 o 

n 30 

n g -n 20 

n g -n 30 

N 2 o N 30 

7051 

-22.57 

-22.08 

-18.49 

-0.48 


-3.59 

7109 

-22.59 

-22.09 

-18.50 

-0.50 


-3.59 

7062 

-32.47 

-32.64 

-30.17 

0.16 


-2.46 

7110 

-30.86 

-31.19 

-28.60 

0.33 


-2.59 

7082 

-12.58 

-12.54 

-8.62 

-0.04 

-3.96 

-3.91 

7084 

-24.90 

-24.46 

-20.69 

-0.44 

-4.21 

-3.77 

7114 

-24.84 

-24.46 

-20.70 

-0.38 

-4.15 

-3.77 

7085 

-29.68 

-29.92 

-26.44 

0.24 

-3.24 

-3.48 

7115 

-29.58 

-30.08 

-26.65 

0.50 

-2.93 

-3.43 

7086 

-20.78 

-19.56 

-17.97 

-1.22 

-2.81 

-1.59 

7885 

-20.83 

-19.56 

-17.97 

-1.27 

-2.86 

-1.59 

7112 

-17.40 

-15.68 

-12.73 

-1.72 

-4.67 

-2.95 

7921 

-29.47 

-27.75 

-25.60 

-1.72 

-3.88 

-2.15 


Units are meters 


(1) ADJUST20 ; 

(2) ADJUST40 : 

(3) ADJUST65 : 

(4) ADJUST71 : 

(5) ADJUST81 ; 


48955 observations 
43271 observations 
48955 observations 
43271 observations 
48955 observations 


(include geophysical data) 
(exclude geophysical data) 
(include geophysical data) 
(exclude geophysical data) 
(include geophysical data) 


The first two solutions were developed based on the same setup as 
far as weighting, removal of high frequency content of the anomalies 
etc. are concerned, as it can be seen from Table 7. The difference 
between them is the exclusion of geophysically predicted anomalies in 
the second set. 


The second pair of solutions differs from the first by the use of the 
improved elevation data (TUG87), and the use of a different weighting 
scheme which will be discussed in Section 5.3.3. Again this pair 
consists of one solution where geophysical data are included in the 
adjustment, and one where they are excluded. 
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Finally, in the last set the systematic reduction <5g h 2 discussed in 
Section 3.1.1, was applied to the anomaly data, in contrast to the 
previous solutions. This solution reflects the optimum procedure for the 
analysis of surface gravity data, suggested by this study. 

The first four normal equation sets already have been used in 
combination solutions performed at NASA/GSFC. 

The majority of the comparions and tests to be presented in the 
sequel refers to solutions (1) and (2). However, as it can be seen from 
Table 7 and as it will be discussed in the following, the differences 
between corresponding solutions (using the same number of 

observations), are rather minor, so that the conclusions drawn from 
these comparisons pertain also to the cases of (3), (4) and (5). 

The evaluation of the quality of our results will be mainly based on 
comparisons with two geopotential models. The GEMT1 [Marsh et al, 
1987] which is a satellite only field complete to degree and order 36, 
and the OSU86F [Rapp and Cruz, 1986b] which is a combined high 
degree geopotential model complete up to 360. 


5.3.1 Comparisons With Existing Geopotential Models 

Before we attempt any comparisons of our results with satellite only 
or combined solutions, it is important to compare the input of our 
adjustment, i.e. the surface gravity data, with the values implied for 
these data from the satellite solutions. Such a comparison will indicate 
the agreement between the gravity field as this is sensed by the 
satellites and the gravity field which is implied by the surface 
measurements. In addition it will provide us with a fairly good idea of 
what we should expect from the comparisons of our results with satellite 
only solutions, since the output of any adjustment primarily reflects the 
characteristics of its input. For that purpose the following comparisons 
were performed: 
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(a) A set of 48955 l*xl* anomaly values were computed as (see Table 6) 

dg in =dg-hf 5 (5.15) 

These values essentially represent the input of our adjustment since 
the systematic correction terms (ellipsoidal corrections etc.) are of 
very small magnitude. 

(b) A global set of l'xl* mean anomalies was formed using the GEMT1 
coefficient set complete from N mln = 2 to N max = 36. These anomalies 
are denoted Ag GEM . 

The differences between Ag GEM and dg tn were then evaluated in the 
areas where the June 1986 field contains an anomaly estimate (Figure 8). 
In Figure 21 the locations where the absolute value of such differences 
is larger than or equal to 10 mgals are shown. From this figure it 
becomes apparent that systematic discrepancies between the satellite 
implied and terrestrial field exist in large continental areas. Since the 
extent of many of these areas exceeds the minimum wavelength 
recoverable by an expansion up to 36, it is expected that the results of 
terrestrial only solutions will reflect these inconsistencies. 

It is also interesting to examine if the above discrepancies are 
related to geophysically predicted data in the June 1986 field. The 
results of a corresponding comparison peformed only over the 43271 
locations of Figure 10 are illustrated in Figure 22. Although some of 
these inconsistencies are associated with geophysically predicted data 
(USSR, China), Figure 22 indicates that such inconsistencies exist also in 
areas where observed gravity data are available (Africa, South America). 

With the above as a background, we next performed the 
corresponding comparisons between Ag GEM and the adjusted values of 
the gravity anomalies obtained from the solution ADJUST20 and 
ADJUST40 (see Table 7), denoted Ag A20 and Ag A40 respectively. The 
results from these comparisons are illustrated in Figures 23 and 24 
respectively. 








Figure 21. Locations of the 17082 l*xl* Blocks Where IAg GEM - dgin' > 
10 mgals (Geophysical Data are Included in the Comparison). 
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Figure 22. Locations of the 14412 l*xl* Blocks Where IAg GEM - dg in l > 
10 mgals. (Geophysical Data are Excluded from the 
Comparison). 
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Figure 23. Locations of the 9745 l*xl* Blocks Where IAg GEM - Ag A20 I > 
10 mgals (Locations Corresponding to Geophysical Data are 
Included in the Comparison). 
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Figure 24. Locations of the 7545 l’xl* Blocks Where *Ag GEM - Ag A40 • > 
10 mgals (Locations Corresponding to Geophysical Data are 
Excluded from the Comparison). 
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There are two main observations to be made with respect to these 
figures. 

(i) 

As expected the inconsistencies between the satellite implied and 
terrestrial fields, that were observed over extended geographical regions 
are also reflected (with remarkable exactness) in the results of 
terrestrial only solution (compare Figures 21 and 23 or 22 and 24). 

(ii) 

The behavior of the least squares adjustment in areas where the 
inconsistencies between the satellite implied and terrestrial fields are 
localized is rather interesting. Typical examples of such cases occur in 
oceanic regions. From the examination of Figures 19, 21 and 23 it 
becomes apparent that the discrepancies between the two sources of 
information in such cases are closely related to the residuals of the 
terrestrial only solution. In fact it appears that the terrestrial 

adjustment yields adjusted anomalies that are closer to the satellite field 
implied anomalies than to the surface observations! 

Before we attempt to explain why such behavior is observed we will 
present a comparison that clarifies the issue. We consider the residuals 
of the solution ADJUST20, denoted v 2 o* and we form the following 
differences: 


a = AjJgem - dg jn 


(5.16) 


and 


b = v 20 - a (5.17) 

Note that (see also Table 6 and equation (5.15)) t 
5 = dgj n + v 2 o - AgcEM ~ Ag A20 - Ag GEM 


(5.18) 
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hence, apart from the small systematic correction terms, b represents 
the differences between the adjusted anomalies from ADJUST20 and the 
GEMT1 implied anomalies. 

Based on the values a and b we then identify the locations where 

lal > 10 mgals and Ibl < 10 mgals (5.19) 

and where 

lal > 10 mgals and Ibl >10 mgals. (5.20) 

Locations where the conditions (5.19) are simultaneously fulfilled indicate 
the cases where, although the input of our adjustment differs from the 
GEMT1 implied anomaly substantially (with respect to the threshold value 
of 10 mgals), the output is closer to the satellite implied anomaly than to 
the terrestrial one (again with respect to the 10 mgal threshold value). 
The cases where the opposite occurs are indicated by the simultaneous 
fulfillment of conditions (5.20). Figures 25 and 26 illustrate the 
geographical distribution of the blocks where conditions (5.19) and (5.20) 
are fulfilled respectively. These figures verify that the speculation 
made in (ii) above about the behavior of the terrestrial solution in the 
cases of isolated (localized) discrepancies between the satellite and 
terrestrial information is indeed a fact. It remains now to explain how 
the terrestrial solution (which has no information about the satellite 
implied field) can yield in these cases results that agree better with the 
satellite information than with its own input. The explanation is rather 
obvious. The least squares adjustment does not yield values agreeing 
better with the satellite implied anomalies; what it does though for these 
isolated blocks, is that it yields values that are in agreement with the 
anomaly values of the broader neighborhood of these blocks (which 
happen to be in good agreement with the satellite implied anomalies). 
The underlying cause of such behavior of the least squares adjustment 
is the sampling interval in conjunction with the maximum estimated 
degree of the expansion. The l*xl* interval provides enough data 
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Figure 25. Locations of the 10149 l*xl* Blocks Where Conditon <5.19} us 
Fulfilled. (Locations Corresponding to Geophysical Data are 
Included in the Comparison). 



Figure 26. Locations of the 6933 l*xl° Blocks Where Conditon (5.20) is 
Fulfilled. (Locations Corresponding to Geophysical Data are 
Included in the Comparison). 
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redundancy to enable the least squares adjustment to isolate (by 
assigning them large residuals), observations which do not comply with 
the more general pattern implied by their broader region (outliers). An 
illustrative example explaining the above can be given if we consider a 
least squares fit of a straight line into sampled data (Figure 27), for 
two different sampling patterns. In Figure 27a all the data sampled are 
used. Point P, which does not comply with information implied by its 
neighboring points will be characterized by a large residual. If now 




(a) (b) 

Figure 27. The Least Squares Adjustment and the Data Sampling. 


the previous observations are averaged over some interval d and the 
same fit is attempted (Figure 27b), the point P which represents the 
mean value in the region of P will be characterized by a much smaller 
residual and the quality of the fit in the region of P will not differ 
much from the quality of the fit on the other points. On the other 
hand if we keep the sampling as in Figure 27a but we start increasing 
the degree of the polynomial to be fitted to the data, the corresponding 
residuals for P will decrease eventually making P appear as a valid 
observation. 

From the above discussion it becomes clear that there is an 
additional good reason for preferring the l*xl* sampling with the 
removal of the high frequency content from the data, over the 
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alternative use of 5*x5* averages. If 5*x5* averages were used we 
would have no way to identify very localized inconsistencies between the 
satellite and terrestrial solutions. 

Up to now all we have determined is the existence of inconsistencies 
between the satellite implied field and the terrestrial one. The next 
step, of course, is to try to determine which of the two sources 
provides a more accurate estimate of the gravity field. To facilitate the 
investigation in this direction we have to employ comparisons with 
external sources of information. We will distinguish between the case of 
the oceanic regions where the inconsistencies found were more of local 
character, and the continental areas where more regional problems were 
experienced. In the oceanic regions, anomalies obtained from altimetric 
observations will provide the external source, while in the continental 
regions (which will be discussed in the next section) Doppler derived 
undulations will be used as independent source of information. 

Figure 28 illustrates the locations of the 6833 l*xl* blocks where 
(see Table 6) 

Idg - Ag alt l > 10 mgals (5.21) 

with &g a it denoting the l*xl* anomaly implied from the altimetric 
observations. These anomalies were estimated from the combined 
GEOS3/SEASAT altimeter data [Rapp, 1985]. In contrast Figure 29 
illustrates the locations of the 3574 l*xl* blocks where (see Table 6) 

lAgcEM — (Ag a it — hf s ) I > 10 mgals (5.22) 

Figure 29 shows a good agreement between the GEMT1 implied anomalies 
and the altimetry derived anomalies. Hence, the discrepancies shown in 
Figure 28 are most probably due to erroneous ship-track measurements. 
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Comparisons as the ones discussed above provide information about 
the geographical distribution and the number of the l°xl* blocks 
where inconsistencies between the satellite-implied and the terrestrial 
data are observed. However, they do not provide information about 
the magnitude of these discrepancies (other than the crude information 
provided by the threshold value used). In order to obtain such 
information we have plotted global undulation and anomaly difference 
maps between the satellite only field (GEMT1) and the terrestrial only 
solutions. 

In Figures 30 and 31 the gravity anomaly and undulation 
differences between GEMT1 and ADJUST20 are illustrated. Note that 
examination of such plots should always take into account the data 
coverage used in the development of the terrestrial only solutions, in 
order to discriminate between areas where large discrepancies occur 
due to lack of terrestrial data (South Pacific, Antarctic region) and 
areas where systematic inconsistencies between satellite and terrestrial 
data exist (Soviet Union, Mongolia, China, Africa and Brazil). 
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Figure 30. Gravity Anomaly Differences: GEMT1 Minus ADJUST20 Based 
on a 5*x5* Grid (Contour Interval is 4 mgals). 
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A comparison in terms of undulations between GEMT1 and ADJUST40 
is shown in Figure 32. As it can be seen from this plot, large 
oscillations of the terrestrial only solution occur in the extended gap 
generated by the exclusion of geophysical data. In addition a 

comparison of Figures 31 and 32 indicates that ADJUST40 is 
characterized by more long wavelength differences with respect to 
GEMT1 than ADJUST20. This should be expected since the exclusion of 
geophysical data, which cover extended regions, degrades the 

estimability of long wavelength features of the gravity field. 

In Figures 33 and 34 the corresponding comparisons are given 
between the OSU86F combined solution (up to N max = 36) and the 
solution ADJUST20. Although these plots are very similar to 30 and 31 
respectively, there is an important difference to be noted from the 
comparison of 30 with 33 (or 31 to 34). The differences between 
ADJUST20 and OSU86F have less power in the high frequencies as 
compared to the differences between ADJUST20 and GEMT1. 
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Figure 33. Gravity Anomaly Differences: OSU86F Minus ADJUST20 Based 
on a 5*x5* Grid (Contour Interval is 4 mgals). 
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LONGITUDE 



Figure 34. Undulation Differences: OSU86F Minus ADJUST20 Based on a 
5*x5* Grid (Contour Interval is 4 m). 


Such behavior should be expected since GEMT1 fails to recover 
high frequency features of the gravity field, while the opposite is true 
for OSU86F and ADJUST20. In Figure 35 this aspect is illustrated in a 
better way. In this figure the percentage differences of PGS3041 (a 
pre GEMT1 satellite only solution), and ADJUST20, with respect to 
OSU86F are plotted. At low degrees (<6) OSU86F and PGS3041 show an 
excellent agreement which is expected since OSU86F has used GEML2' 
[Lerch et al, 1982a] to control these coefficients, and PGS3041 is not 
expected to have large discrepancies with GEML2'. 

For the same degrees, in the case of ADJUST20 the percentage 
differences are larger, a fact which expresses the inability of the 
terrestrial only solution to recover long wavelength features of the 
geopotential. However, at degree about 20, the picture changes, and 
PGS3041 departs from OSU86F, while the discrepancies of ADJUST20 
with OSU86F start to drop. This is due to the inadequacy of a 
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Figure 35. Percentage Differences Between OSU86F Versus ADJUST20 
and OSU86F Versus PGS3041. 

satellite solution to recover high frequency features of the 
geopotential. Figure 35 can be used to describe in a compact and 
clear manner the whole idea and purpose behind the combination 
solutions. It also suggests that when Figures such as 30 or 31 are 

examined, very long wavelength biases should not be considered a 
problem, since they merely reflect the inability of the terrestrial only 
solutions to recover low degree coefficients. However, systematic 
differences such as the ones occuring in Africa, constitute a problem, 
since their extent corresponds to wavelengths which should be 
recoverable by terrestrial only solutions. 
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5.3.2 Systematic Inconsistencies Between Satellite Only and Terrestrial 
Only Solutions 

One of the continental regions where systematic inconsistencies 
between the satellite and terrestrial solutions are observed is Africa. 
We have decided to further investigate the problem in this area for 
two reasons: 

(a) The discrepancies in this area are quite large (e.g. 24 m undulation 
differences are observed in Figure 31), and 

(b) A fairly well distributed set of Doppler Stations covering the area 
in question was available, providing an independent standard for 
comparison. 

Large systematic differences occur also in Asia (see Figure 31), 
however no redundant information for comparisons is available in these 
regions. 

In Figure 36 a more detailed anomaly difference contour map is 
given for the area in question. It can be seen from this figure that 
the differences reach values of 24 mgal. 

On the other hand, Figure 37 shows that the undulation differences 
in this area reach values of -24 m. Figure 37 shows also the locations 
of some of the 32 Doppler stations that were used in the following 
comparison: 

(a) We started with a set of 32 "Doppler undulations" denoted by N 0 
and evaluated by 

N 0 = h 0 - H (5.23) 

where h 0 (the Doppler derived ellipsoidal height) was first 
converted from the Doppler reference frame to a geocentric, 
properly scaled reference frame using the translation and scale 
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Figure 36. Gravity Anomaly Differences: GEMT1 Minus ADJUST20 Based 
on a 2*x2* Grid (Contour Interval 2 mgals). 
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Figure 37. Undulation Differences: GEMT1 Minus ADJUST20 Based on a 
2*x2* Grid (Contour Interval 4 m). Superimposed are 

Selected Doppler Station Locations. 
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Table 14. Undulation Comparisons at Doppler Stations in Africa. 


SELECTED DOPPLER STATIONS IN AFRICA 


N* » DOPPLER UNDUL. - UNDUL. IMPLIED BY OSU86F (NMIN-37.NMAX-360) 
NT1 l UNDUL. IMPLIED BY GEMT1 (NMIN=*2,NMAX»36) 

NA20 l UNDUL. IMPLIED BY ADJUST20 (NMIN-2.NMAX-36) 

DN1 - N* - NT1 
DN2 - N* - NA20 

DN3 - DOPPLER UNDUL. - UNDUL. IMPLIED BY OSU86F (NMIN-2.NMAX-360) 


STN. 

LAT 

LON 

H{ ELL) 

N* 

NT1 

NA20 

DN1 DN2 

DN3 

305 

-4.86634 

29.63628 

907.830 

-13.77 

-12.62 

-1.66 

-1.16 -12.11 

0.03 

333 

-1.16821 

36.94086 

1518.052 

-16.49 

-20.35 

-5.04 

3.87 -11.45 

0.99 

334 

-1.13384 

31.38364 

1314.780 

-13.84 

-13.68 

2.53 

-0.17 -16.38 

-1.13 

345 

-0.57266 

34.85817 

1999.122 

-13.58 

-17 * 16 

0.34 

3.57 -13.92 

0.87 

358 

0.51367 

25.18820 

392.320 

-17.16 

-17.52 

-0.48 

0.35 -16.68 

-2.43 

370 

2.26306 

37.90152 

1406.962 

-15.37 

-19.55 

-1.77 

4.18 -13.60 

0.35 

373 

2.71574 

36.68937 

368.812 

-13.37 

-17.38 

2.59 

4.01 -15.95 

0.36 

375 

3.11495 

35.60518 

486.552 

-11.68 

-15.76 

5.50 

4.08 -17.17 

1.15 

379 

3.51501 

38.55354 

678.762 

-16.49 

-18.57 

-0.53 

2.09 -15.95 

-2.15 

381 

3.84124 

11.52424 

767.720 

15.67 

14.03 

23.33 

1.64 -7.67 

0.36 

394 

4.59185 

13.68789 

715.800 

13.32 

9.94 

22.39 

3.38 -9.08 

0.93 

395 

4.75553 

31.63069 

630.020 

-8.32 

-11.17 

12.16 

2.85 -20.48 

1.90 

440 

7.17102 

38.30381 

2202.210 

-9.15 

-11.45 

8.78 

2.30 -17.93 

-1.54 

450 

7.69683 

28.00029 

445.630 

-0.03 

-7.30 

19.58 

7.28 -19.60 

2.40 

462 

8.72038 

38.88549 

1904.012 

-9.67 

-9.53 

9.07 

-0.14 -18.74 

-3.02 

468 

8.98466 

38.78894 

2316.122 

-6.19 

-8.94 

9.54 

2.75 -15.73 

0.23 

471 

9.08775 

36.58082 

2159.125 

-2.48 

-6.37 

14.19 

3.88 -16.67 

3.26 

510 

10.99886 

27.83403 

429.510 

0.99 

-2.47 

22.76 

3.47 -21.77 

0.11 

512 

11.08385 

39.71495 

1853.680 

-8.09 

-7.14 

7.73 

-0.95 -15.82 

-1.34 

530 

11.93986 

41.45189 

409.380 

-9.44 

-9.04 

1.86 

-0.39 -11.30 

0,14 

535 

12.51986 

37.43329 

1959.740 

-1.99 

-2.02 

14.61 

0.03 -16.61 

1.29 

538 

12.73172 

43.11202 

53.260 

-9.75 

-10.89 

-3.81 

1.14 -5.95 

2.73 

541 

13.10689 

33.85169 

428.375 

-1.68 

-0.55 

19.57 

-1.13 -21.24 

-0.18 

543 

13.43213 

22.43807 

815.550 

6.67 

8.83 

24.81 

-2.16 -18.14 

-2.67 

547 

13.56639 

25.22098 

802.260 

8.04 

6.08 

25.91 

1.96 -17.87 

0.38 

575 

15.14766 

36.10997 

481.653 

2.27 

1.61 

18.83 

0.66 -16.55 

1.31 

584 

15.53997 

32.48451 

397.589 

2.96 

2.92 

23.29 

0.04 -20.33 

-0.55 

585 

15.61135 

32.53878 

404.400 

2.57 

3.02 

23.33 

-0.45 -20.76 

-1.05 

594 

16.16405 

32.75606 

430.011 

3.85 

3.85 

23.71 

0.00 -19.86 

-0.64 

622 

17.93626 

30.96552 

304.768 

8.27 

8.31 

25 • 75 

-0.03 -17.47 

0.75 

661 

19.50155 

37.26091 

4.417 

3.77 

3.46 

15.84 

0.30 -12.07 

1.36 

662 

19.58522 

33.31354 

362.118 

7.46 

7.83 

23.19 

-0.36 -15.72 

0.64 


parameters derived by Boucher and Altamimi (1985) 

(b) We removed the high frequency contribution to these undulations 
using the OSU86F model [Rapp and Cruz, 1986b] from N m | n = 37 to 
N max = 360 

(c) We computed the undulations implied by GEMT1 and ADJUST20 for 
the locations of these stations. 

The results from these computations are given on Table 14. The last 

column of this table shows the differences between N 0 and the 
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undulations implied by OSU86P (from N roin = 2 to N max = 360), as an 
indicator of the agreement of Nq with undulations implied by a high 
degree geopotential model. 

As it can be seen from this table, the undulations implied by 
GEMT1 are in good agreement with the Doppler derived undulations 
(the RMS of the differences Np - N GEM ti is 2.6 m), while the 

undulations implied by ADJUST20 are systematically biased with respect 
to the Doppler undulations (the RMS difference Np - N A0JUST2 o is about 
16 m). 

This comparison indicates that the problem in Africa is due to 
erroneous surface gravity values. This was also verified by another 
comparison performed between 10* equal area mean gravity anomalies 
obtained from Satellite-to-Satellite Tracking (SST) techniques given by 
Kahn et al (1982) with the corresponding 10* equal area mean values 
computed from the June 1986 l*xl* anomalies by area averaging. The 
results of this comparison are given on Table 15. 


Table 15. 10* Equal Area Gravity Anomaly Comparisons in Africa. 


Border of 
10* ea. Block 

10* e^. Meai 
Ag (mga] 

Anomaly 

LfQ 

El 

El 

^lU 


SST 

GEMT-l 

JUNE86 

B9 

-15* 

5* 

10:8825 
31 : 1925 

21: 0375 
41:3475 

10.9 

7.1 

8.2 

7.2 

-8.9 

22.7 


The 10* e.a. anomalies implied by GEMT1 (computed by area 
averaging of the l*xl* equi-angular implied anomalies) agree well with 
the SST derived anomalies, while the terrestrial anomalies show large 
differences. In both the cases of June 1986 as well as of GEMT1, the 
l*xl* equi-angular blocks used to define the 10* equal area gravity 
anomaly were determined by rounding the 10* block longitude borders 
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to the nearest degree. 

The regional inconsistencies between the anomalies implied from 
satellite only geopotential models and the terrestrial measurements were 
identified also in previous investigations [Kahn et al, 1982], [Mainville 
and Rapp, 1985]. Possible causes for these biases (in the cases where 
observed gravity data exist) are: 

(i) Systematic biases in gravity base stations 

(ii) Systematic errors introduced in the collection or compilation of 
the gravity data (gravity formula etc) 

(iii) Systematic biases of the vertical datums of these regions (see 
also Section 2.3.5). 

In the region of Africa, examination of the differences between the 
TUG86 and the TUG87 l # xl* mean elevations has not shown any 
systematic differences in the region that could account for the anomaly 
differences detected. However, this is only to be interpreted as an 
indication since we have no information about the elevations used for 
the compilation of the l # xl* mean anomalies in that region. 

On the other hand, in the areas of geophysically predicted 
anomalies, the biases may be due to inappropriate modeling and 
estimation techniques used to estimate these values. 

Finally, it should be mentioned that the regional inconsistencies 
detected here have also been experienced in previous combination 
solutions. Figure 26 is very similar (over continental regions) with the 
Figure 11 given in [Rapp and Cruz, 1986a, p. 46] which illustrates the 
locations of residuals exceeding 7 mgals in the OSU86D combined 
solution. 
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5.3.3 Downweighting of Selected Terrestrial Gravity Data 

Three alternative procedures can be used in order to compensate 
for the regional inconsistencies discussed in the previous section. 

(1) 

In areas where it is identified that the problem is caused due to 
vertical datum inconsistencies, introduce additional parameters in the 
modeling as it was discussed in Section 2.3.5. 

(2) 

In areas where redundant information is available (e.g. Doppler 
undulations), introduce additional observations in the adjustment. 

(3) 

Downweight the anomalies in regions where it is identified that the 
satellite information provides a more accurate description of the 
gravity field. 

Method (1) suffers the disadvantage that the existing correlations in 
the normal matrix may cause contamination of other parameters instead 
of resolving vertical datum biases. In any case the use of such a 
method requires thorough investigation that is beyond the scope of 
this study. 

Method (2) is limited by the fact that redundant information is 
seidomly available. For example, no Doppler undulations are available 
in the regions of the Soviet Union or China where large regional biases 
were detected. In addition, introduction of sparse point undulation 
data in a global adjustment has to be studied carefully. 
Parameterization of these observations should not include the very low 
degree harmonics since the sparse distribution of these data will 
degrade the estimability of long wavelength features of the gravity 
field. In an experimental solution where a limited set of point Doppler 
undulations was used in Africa as additional observations, a reduction 
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of the undulation differences (GEMT1 minus terrestrial solution), by 4 
m was obtained. However, in this solution only coefficients above 
degree 20 have been parameterized in the undulation observations and 
a very Bmall relative weight (0.01) was used in the combination of the 
normal equations obtained from the undulations with the normals 
obtained from the anomalies. 

Such a solution was made only for experimental purposes since we 
recognize that the use of heterogeneous data (both in terms of nature 
and in terms of distribution), is a problem that could not be addressed 
in the time frame of this work. 

Consequently, we have decided to attempt to compensate for the 
inconsistencies using method (3). The regions where the data would be 
downweighted were defined through the following comparisons: 


(a) 

| Ag 1 * x i«osue6F 3 ^ 36 “ 

Agi» x i* 

ADJUST *0 1 > 

10 mgals 

(5.24) 

(b) 

1 Agj* x i *0SU86E 2 ^ 36 - 

Agi'xi* 

ADJUST 40 1 > 

10 mgals 

(5.25) 

The 

locations of these 

blocks 

are shown 

in Figures 

38 and 


respectively. 



Figure 38. Locations of the 2550 l*xl* Blocks Where I Ag x » xX ‘osuecF 
~ Ag x » x j *aojust 2 o I * 10 mgals. 
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Figure 39. Locations of the 1985 l*xl* Blocks Where I Agi * x i •osue«£ 2 -» 3 « 
- Ag 1 » xl * A0JUST4O l > 10 mgals. 


The OSU86E/F implied l*xl* anomalies used in the comparisons 
contained the part of the spectrum only up to N max = 36. 

Two files containing "flags" indicating these locations have been 
constructed (see also Section 5.2.1). The standard deviations of the 
anomalies in these locations were defined by 




max(30, 2.5<r, j ( JUN8 «) ) 
min (60, 2.5<r, j ( JUM86 ) ) . 


(5.26) 


In the case corresponding to Figure 38, the area average standard 
deviation of the 46405 blocks which do not fulfill (5.24) (these standard 
deviations are assigned according to (5.11)), is 28.8 mgals. The area 
average of the standard deviations of the 2550 blocks (assigned 
according to (5.26)) is 44.2 mgals. Hence, a downweighting of about 2.5 
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is performed for the areas where inconsistencies occur. 

Two solutions performed according to this weighting scheme are 
the ADJUST65 (including geophysical data) and ADJUST71 (excluding 
geophysical data) (see also Table 7). The comparison of these solutions 
to the corresponding solutions ADJUST20 and ADJUST40 gave the 
overall results shown on Table 16. 

As it is expected the downweighting affects more the solution 
where geophysical data are excluded because in that case the total 
number of available observations is substantially smaller than when 
geophysical data are included. 

In addition, the comparison of ADJUST65 against the GEMT1 field in 
terms of global undulation and anomaly maps did not show substantial 
changes over the corresponding comparison of ADJUST20. This is due 
to the fact that the effect of downweighting cannot become apparent 
unless redundant information is introduced for the regions in question 
in the combination solution with the satellite normal equations. 


Table 16. Overall Comparison Between Terrestrial-Only Solutions 
Without and With Downweighting of Selected Anomaly Data 



ADJUST20 

ADJUST40 


versus 

versus 


ADJUST65 

ADJUST71 

Average % Diff 

8.67 

10.49 

RMS Undulation Diff (m) 

0.60 

3.67 

RMS Anomaly Diff (mgals) 

1.16 

6.29 


We recognize here that the procedure adopted to compensate for 
the regional inconsistencies between satellite-implied and terrestrial 
anomaly data, is not an optimum one. These problems should be 
considered in a more systematic manner in future updates of global 
anomaly data bases. 
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5.3.4 Computation of Anomaly Degree Variances 

Concluding this chapter an important issue related to the 
computation of anomaly degree variances from coefficients obtained 
from a least squares adjustment of an incomplete set of gravity 
anomaly data will be discussed. 



Figure 40. Anomaly Degree Variances Implied by GEMT1 and ADJUST20. 


DEGREE 


Figure 41. Anomaly Degree Variances Implied by GEMT1 and ADJUST40. 


The anomaly degree variances implied by a given potential coefficient 
set were computed according to 
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c n = y a (n-l) 2 Z (c£ + §L) (5.27) 

m-o 

where 



The values of GM, a used were the ones given in (3.3), while the even 
degree zonal coefficients of the normal potential (J 2 , J 4 , J a ), used to 
define C* m (see equation (3.10)), were the ones given in (3.13). 

Figures 40 and 41 show the anomaly degree variances implied by 
the satellite only field GEMT1 and the terrestrial only fields ADJUST20 
and ADJUST40. 


Examining Figure 40 we observe that the solution ADJUST20 gives 
anomaly degree variances which are in good agreement with those 
implied by GEMT1. Also, as it should be expected the terrestrial 
solution appears to have more power at higher degrees (above degree 
about 20), than the satellite one. 

On the other hand, a first look of Figure 41 suggests that 
ADJUST40 yields coefficients totally unrealistic. However, the awkward 
shape of the spectrum implied by ADJUST40 is due to the fact that in 
the computation of degree variances, the assumption is inherent that 
the potential coefficients are uncorrelated. In the case of ADJUST20 
the correlations between the coefficients which are neglected are small 
as it can be deduced from Figure 15a or 15c. However in case of 
ADJUST40 as it can be deduced from Figures 15b or 15d, the 
correlations between the coefficients cannot be omitted. 


Degree variances, in theory, cannot be defined for any set of 
coefficients estimated from gridded data (point or mean values) of 
finite gridsize, even if the data grid is complete, according to what 
was discussed in Sections 4.2.1 and 4.2.2. However as it can be seen 
from the above the omission of the existing correlations between the 
coefficients effects the results severely only in cases where the 
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coverage is quite limited. 

Computation of "degree variances" can be done if we diagonalize 
the normal matrix associated with the coefficients. However, the base 
functions that will be defined from the diagonalization will be linear 
combinations of surface harmonics not necessarily of the same degree . 
Hence, we would be able to compute the contribution to the anomaly 
from different coefficients but these coefficients will not necessarily be 
of the same degree. 


CHAPTER VI 


SUMMARY - CONCLUSIONS - RECOMMENDATIONS 


The estimation of harmonic coefficients of the geopotential from 
surface gravity data was reexamined. The ultimate purpose was to 
investigate optimum modeling procedures and estimation techniques, in 
order to incorporate the information implied by the most up to date 
terrestrial gravity data, in combination solutions currently developed at 
NASA/Goddard Space Flight Center for use in the TOPEX/POSEIDON 
mission. 

The modeling adopted is based on a form of the fundamental 
boundary condition of the linerarized Molodensky boundary value 
problem, which accounts for the effects of the earth’s ellipticity to the 
order of the eccentricity squared. 

The surface anomaly data were corrected for the effects of the 
atmosphere and the earth’s ellipticity. The ellipsoidal correction terms 
were computed from an existing geopotential model, and applied to the 
surface anomaly data prior to the adjustment. In addition, examination 
of the procedures used for the numerical realization of the Molodensky 
surface free-air anomalies, lead to the conclusion that the term 
corresponding to the second order vertical gradient of the normal 
gravity (<5g h 2) is neglected in the evaluation of the surface anomalies. 
This term causes long wavelength systematic errors in the undulations 
that have an RMS magnitude of about 22 cm, with maximum values 
reaching 1.80 meters. 
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Consequently, the additional correction (fig h 2 ) was applied to the 
surface anomalies to compensate for this effect. 

The observation equations used to relate the unknown harmonic 
coefficients to the gravity data (properly reduced for the above 
systematic effects), consider the data lying on the physical surface of 
the earth. Hence, analytical continuation of the surface data to a 
bounding surface, is not required for the solution. 

The estimation of harmonic coefficients from the surface anomalies 
was peformed using a least squares minimum variance estimator as 
proposed by Rapp (1967). The full normal matrix was formed and 
inverted to obtain the least squares solution. The advantage of this 
technique is that it does not require a complete set of anomaly data, 
and that it does not involve any assumptions regarding the 
orthogonality of area-mean samples of surface spherical harmonics. 
Hence, in contrast to the quadrature formulas, it enables exact 
recovery of a set of harmonic coefficients from simulated (noiseless) 
anomaly data computed from these coefficients. 

An analytical comparison between the quadrature formulas and this 
technique, resulted in an alternative way which is proposed here for 
the definition and numerical evaluation of the quadrature formulas 
given in equation (4.46). Such a definition does not require the use of 
optimal quadrature weights. 

In the application of the least squares adjustment procedure, an 
alternative approach to compensate for aliasing effects was proposed 
and used. Namely, instead of predicting mean anomalies of gridsize 
that corresponds to the maximum degree and order of the coefficients 
being estimated (AX = rr/N max ), we use the original l*xl* anomaly data 
and remove from them the frequency content beyond N max . The high 
frequency content of the data is computed from an existing high 
degree geopotential model. 
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This procedure avoids the prediction of mean anomalies from the 
l*xl* data, which would introduce additional approximations, and 
provides a better control over the frequency content of the data. In 
addition, the use of more dense data in low degree expansions enables 
the distinction of isolated erroneous data values, which are 
characterized by large residuals. From this point of view, such low 
degree expansions can be used as efficient tools for data editing in the 
compilation of global anomaly fields. 

The comparison of low degree (N max = 36) global gravity models 
obtained in this study from terrestrial data alone, with the satellite 
derived models, has indicated extended geographical regions (USSR, 
China, Africa) where systematic inconsistencies between the two sources 
of information occur. In the case of Africa comparisons with Doppler 
derived undulations verified that the problem lies with the terrestrial 
data. In an attempt to restrict the influence of the data from such 
regions in the combined solution, a downweighting of these anomalies 
was performed. It is recognized that this procedure is not really an 
optimum one and in future compilations of global anomaly fields special 
emphasis should be given in the improvement of the gravity material 
for these areas. Also, the techniques used for the estimation of 
gravity anomalies through geophysical correlation should be reexamined 
so that more reliable gravity anomaly estimates can be obtained from 
the application of these procedures. 

The alternative use of gravity anomalies reduced to the ellipsoid 
through the application of terrain corrections was tested on a limited 
area. This procedure did not yield satisfactory results. Additional 
investigation is required in order to clarify the use of terrain 
corrections in the development of global geopotential models. 

The nature of the topic of this study has introduced us to a 
number of related problems for which detailed investigation could not 
be conducted, inside the limits of this work. The estimation of error 
properties for global anomaly fields is certainly one of the problems 
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that should be investigated in the future. The study by Weber and 
Wenzel (1982) can serve as a starting point for such investigation. 

The implementation of heterogeneous data (gravity anomalies, sea 
surface heights etc) in global geopotential modeling is another topic 
that deserves additional investigation. The ability to incorporate 
different and repeat data in the least squares adjustment, makes this 
estimation procedure a natural selection for solutions to the 
altimetry-gravimetry boundary value problem in mixed domains. 
Solutions as the one performed by Wenzel (1985), should be further 
investigated in an attempt to avoid some of the approximations used 
(e.g. formation of diagonal normal matrix). 

Finally, the usefulness of the modern computational resources in 
the investigations related to global gravity modeling should be 
emphasized. This study serves as a perfect example of an 
investigation that was highly facilitated through the use of a 
supercomputer. Solutions which were carried out in a routine basis 
for this study were avoided 10 years ago due to their computational 
load. It is our sincere hope that these computational resources will 
make the term "computational limitations" to appear obsolete in a few 
years. 



APPENDIX A 


SYSTEMATIC EFFECTS 


In this section the systematic effects arising from three causes are 
examined in terms of numerical values and spatial behavior. These are: 

(a) « h : represents the effect on the boundary condition of the linear 

(first order) term in the Taylor series expansion of 3/3h 
around its "sperhical" value 3/3r, and thus has its origin in 
the Taylor series expansion of the geometrical properties of the 
ellipsoid around their spherical values. 

(b) Eyl represents the effect on the boundary condition of the linear 

term of the Taylor series expansion of the normal gravity y 
and its gradient 3y/3h around their spherical values. Hence, it 
originates from the Taylor series expansion of the dynamical 
properties of the ellipsoid around their spherical values. 

(c) e r : represents the effect on the area-mean gravity anomaly caused 

by neglecting the variation of the geocentric radius over the 
block (omission of the second term in equation (2.99)). 

The analytical expressions for the first two terms are 

e h = e 2 sin0cos0 [—j] (A.l) 
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c 7 = [ 6 J 2 fr P 2 (cos0) - ^~r- sin 2 ©] T (A. 2) 

and the formulas for the computation of area-mean values IE h and IEy 
of these terms are given in (2.111) and (2.112). The corresponding 
formula for * r is derived next to the 0(e 2 ). We have (see equation 
(2.99)) 

!•(♦) = fjj - \ N,(M, + H,j)(* - *,) + ... (A. 3) 


so that 


1 _ _i_ 

«*(♦) r,j 


l 1 + f 


1 e 2 sin2$< r , <•« 




N,(M f 


+ H < j ) (♦ 


*i> +...] 


(A. 4) 


and 




+ H, j)(* - ♦,) + ...] 
(A. 5) 


or 


twf ■ wo + ««i 

with obvious notation. 

However, we have 


Ag(r,0,X) = 


77 l (n— i> (j) ” I (C* m cosmX + S nm sinrnX)P nm (cos0) 

r n = 2 m=o 

(A. 7) 


Setting 


D nm = C* m cosmA + S nm sinmX 


(A. 8 ) 



and 
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Pnm = Pnm( cos ®) 

equation (A. 7) can be written as 




D n0l P 


nm 


or according to (A. 6) 


D 4* 2 

Ag * “7 E (n-l)LM [1 + X«0] £ D n(n P 

a n=2 vr 1j-' m=o 


or 


« g °f2 i (n-DkS-p E 

a n = 2 m = 




GM 

a 2 



DnmX(*)P„ m 


so that 

Ag = Ag* + c r 

where Ag* is the anomaly computed with the geocentric 
referring to the center of the block, and 

c r = e 2 p in f 1 N,(M, + H 1j )GM £ (n-l)(n+2)M-)" £ D nro (* - 

^ r ij n=2 'Mj' m=o 


or 


£ (n-1) (n+2) MM " £ D n-I (* - »,)P r 

n=2 '■ r 1 J J m=o 


(A. 9) 
(A. 10) 

(A. 11) 

(A. 12) 

(A. 13) 
radius 

*i)Pnm 

(A. 14) 
(A. 15) 


For small blocksizes (e.g. l*xl*) we may write 
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(♦ - $t) rad 8 sin($ — = sin«frcos? i - cos$sin$j (A. 16) 

so that (A. 15) becomes 

e r = c cos?, 1 (n-l)(n+2)M-] n Z D nn) sin*P nm 

n = 2 j ' m— o 

- c sin$( Z (n-l)(n+2) U 2 -]" Z D nm cos4»P nm (A. 17) 

n=2 j ' m=o 

However, (see [Rapp, 1984, p. 26, equation 3.78]) 

♦ = Y' + msin2? + ... (A. 18) 

or 

♦ 8 V + |- sin24> + ... (A. 19) 

where V denotes the geocentric latitude = 90* - 9). The last term in 
(A.19) will introduce terms of 0(e 4 ) when substituted in (A.17), hence it 
can be omitted. Accordingly (A.17) becomes 

c r = c cos? j Z (n-l)(n+2)(^-)" Z D nm cos0P nm 

n — 2 vr i j J m=o 

- c sin?, z (n-1) (n+2) MM Z D nm sin0P nm (A. 20) 

n=2 ' r 1 j ' m=o 

From the recursive relations [Moritz, 1980, equation 39-45] 

COS0P nm (cos0) - 2 n+ i Pn+I,m + 2 n +l (A. 21) 

and [Ilk, 1983, equation Z.1.28] 

sin0P nm (cos0) = 2^1 (Pn+I,m+1 ~ P n -i,„,+i) (A. 22) 


we have 
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COS0P nm (coS0) x nm^n+i,m + ynm^n—ijin (A*23) 

sin0P nm (cOS0) - z nm^ n+i,m+i ” w nmPfi— l,m+i (A, 24) 

where 


/ (n-m+1) (n+m+l7 f x 

Xnm - J (2n+l) (2n+3) (a) 

y " m =/^» (n > (b) 

/~LP.±5.t 1 JI n+in+2) ( . 

z nm •/ (2n+l) (2n+3) ( ' C; 

Wnro >/ (2n+l) (2n-l) L) W 


(A. 25) 


Hence, 


®r 


c { cos$,[J a (n-l)(n+2)(^-) n J o D nm x nra P n+1 >( 


+ I (n-l)(n+2)M-)" "z D nm y nm P n — j m ] 

n- 2 Kr ij J m=o J 



2 


(n 1) (n+2) (p ) Dnro z nin^n+i , tn+i 


- z (n-1) (n+2) (p^“) ” Y D nm w nm P n -x, m+l ] } 

m=o lr lj J m=o ’ 1 > 


(A. 26) 


Evaluation of point values of e r can be made using (A. 26), truncated to 
some finite degree N max , from a given set of potential coefficients, 
while for area-mean values (IE r ), the above equation has to be modified 
in exactly the same manner as it was done in the case of e h and ty 
(see equations (2.111) and (2.112)). As in the case of e h and ty , the 
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evaluation of s r (or IE r ) may be made on the surface of the reference 
ellipsoid rather than on the physical surface of the earth without 
significant loss of accuracy. 

For the numerical application , l*xl* mean values of the terms IE h , 
IEy and IE r have been computed from the OSU86F potential coefficient 
set [Rapp and Cruz, 1986b] complete to N max = 180, using the TOPEX 
constants. 

A harmonic analysis was then performed to each one of the global 
sets containing the anomaly correction terms IE h , IEy and IE r . From 
the corresponding harmonic coefficient sets the implied undulation 
corrections were then computed. In Figures 42, 43 and 44 these 
undulation corrections are illustrated. These plots represent point 
values of the undulation on a 5*x5* grid, computed from harmonic 
coefficients complete to N max = 36. 


LONGITUDE 



Figure 42. Undulation Corrections Implied by IE h . (Contour Interval is 
5 cm). 
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With respect to these plots the following comments are to be made 

(a) The term e h is proportional to the meridional component of the 
deflection of the vertical, since 

e h * e 2 sin0cos0-yf (A.27) 

as it can be seen from (A.1). It is thus expected to take large 
values wherever the meridional slope of the geoid is large. An 
example of this behavior is the high observed in Figure 42 in the 
area of the Indian Ocean. 

(b) The term tj is proportional to the undulation as it can be seen 
from (A.2). This explains the similarity of Figure 43 with global 
plots of long wavelength geoids (compare, for example, with 
[Heiskanen and Moritz, 1967, Figure 3-21]). 

Both Figures 42 and 43 are in perfect agreement with the 
corresponding Figures 1 and 2 in [Cruz, 1986, pp. 15, 16] verifying the 
equivalence of the two methods for the computation of the ellipsoidal 
effects. 

(c) Figure 44 indicates that the undulation corrections implied by t r 
are almost exclusively a function of latitude, which is of course 
expected. It also indicates that the effect of t r can be safely 
neglected since the maximum undulation correction implied by this 
term is (in absolute value) about 3 cm. 
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APPENDIX B 


COMBINATION SOLUTION ADJUSTMENT MODELS 


Notation 


X a : theoretical values of parameters 
L a : theoretical values of observables 
L b : observed values of observables 

A g 

L : adjusted values of observables 
X°: approximate values of parameters 
X : adjusted values of parameters 
v : residuals 

L x : a set of observed parameters 

: alteration vector (solution of normal equations) 


(B.l) 


a = X° + X (B.2) 

s = L b + v (B.3) 

Identify: 

L x : potential coefficients obtained from satellite solution (B.4) 
P x : weight matrix associated with L x (B.5) 

L b : area-mean gravity anomaly estimates, dg 1 j , properly 
corrected for systematic errors 
Pf: weight matrix associated with L b 
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(B. 6) 
(B. 7) 


(B. 8) 
(B.9) 
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Assume : 
Set: X° 


I* x 3nd 


uncorrelated 


= L„ 


Case I : Terrestrial estimates of potential coefficients obtained from 
quadrature formulas 

Mathematical Model: 


C 


a 

nm 


1 

47ry (n— 1) 


N-l 

I IPn, 

i -O 


2H—1 

l 


ic 

.IS 


j 

m 

j 

m 


dSi j 


(B. 10) 


or 


X a = G(L a ) 


(B. 11) 


Combination solution can be viewed as the sequential adjustment 


X a = G(L a ) 
L x = X a 


(B. 12) 


which reduces to the case of "observations in two groups with the 
same unknown parameters" [Uotila, 1986, Section 8], as long as we 
identify 


Fi(L?,X a ) = 0 <=> X a - G(L a ) = 0 
F a (Lj,X a ) = 0 <=> L x - X a = 0 

with L, = Lj; L a = L x 

Then, 


(B. 13) 


(B. 14) 


A x 


▼ a 
u 


X a 


_ 

= Lt 


= x v 


= I 


3F, 

3X a 


(B.15) 
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W x = F x (L?,X°) = L x - G(Ij) := W 



B a 




I 


W 2 = F 2 (L 2 ,X°) = L x - L x = 0 


Accordingly, 



M = BP-B T = ft' °J = 

so that 


Mx 1 0 


(B*P*‘Bi) 

0 

.0 Ml 1 . 


0 

Px- 


BiP/Bj 0 


—i 


(B. 16) 

(B. 17) 
(B. 18) 

(B. 19) 
(B.20) 


(B.21) 


(B.22) 


(B. 23) 



yields 


X = -[(B*P* X bI) 1 + P x ] ‘(BiP^Bj) X W 


Also 


v = P X B T K 

where K are the Lagrange multipliers. Hence 


C3- 

Pi* o 


Bi 0 

[Kx] => Vi = Pi l BiKx 

.0 p; 1 . 


.0 I. 

KaJ v x = P x *K a 


Now 


[y = <“ + »> 



yields 

Kx = -(BiPrBj) _1 (X + W) 

K 2 = P x X 

From (B.26) and (B.27) we have 
Vi = -Pi 1 Bi(BiP7 1 Bj)" 1 (X + W) 

v x = X 

Equation (B.29) in view of (B.24) yields 


(B.24) 


(B.25) 


(B.26) 


(B.27) 


(B.28) 

(B.29) 
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v x = -[(BiP/Bj) * + P x ] '(BiPi'Bl) *W (B.30) 

Equation (B.30) is identical with equation (74) in [Rapp, 1984], while 
(B.28) is identical to [ibid, eq. 69] since B x = I in Rapp’s notation. 

Case II : Terrestrial estimates of potential coefficients obtained from 
least squares adjustment. 

Mathematical Model 


dgij 


1 GM 

Affj r^ 


Nmax 

Z 

n = o 




m-o <rt-o 


(B.31) 


The combination solution can be viewed as the sequential adjustment 


L a = F(X a ) 
L x = X“ 


(B.32) 


which again falls into the same category as before [Uotila, 1986, Section 
8] as long as we identify 


Fi (L a , X a ) = 0 <=> L a - F(X a ) = 0 
F 2 (L a ,X a ) = 0 <=> L x - X a = 0 


(B. 33) 


Following the same procedure as in Case I we now have 



B x = I 

Wx = Fx(Lx,X ) = Ly - F(L X ) := W 
A a = -I 


(B.34) 


(B. 35) 
(B. 36) 


(B.37) 
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(B.38) 

(B.39) 



(B.42) 

identical to equations (117) and 
388] once the proper changes of 
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